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Numerical prediction of steady and unsteady tip 
vortex cavitation on hydrofoils 


Paweł Flaszynski, Ph.D. 

Jan A. Szantyr, Prof. 

Krzysztof Tesch, Ph.D. 

Gdansk University of Technology, Poland 


ABSTRACT 


The article presents the numerical method for prediction of tip vortex cavitation generated on hydrofoils. 
This method has been developed in the course of numerical and experimental research described in earlier 
publications. The objective of the research was to design the optimum discrete grid structure for this specific 
computational task and to select the best turbulence model for such an application The article includes 
a short description of the method and a computational example demonstrating its performance. In this 
example the results of numerical prediction of the cavitating tip vortex obtained from two commercial CFD 
codes are compared with experimental photographs taken in the cavitation tunnel in the corresponding 
flow conditions. Altogether nine different flow conditions are tested and analyzed, but only selected results 
are included. The accuracy of the numerical predictions is discussed and the reasons for minor existing 
discrepancies are identified. The unsteady tip vortex calculations are also presented, showing the fluctuations 
of the transverse velocity components predicted for three cross-sections of the cavitating vortex kernel. 


Key words: rotary hydraulic machinery; vortex cavitation; numerical methods 


INTRODUCTION 


The computational method used in this article has 
been developed in the course of detailed experimental and 
numerical research described in [2, 3, 4]. This research was 
at first devoted to numerical prediction of the tip vortices 
generated behind the tips of hydrofoils without cavitation. 
At that stage the objective was to reproduce as accurately 
as possible the complicated velocity field generated by the 
non-cavitating tip vortices. This velocity field is dominated 
by the strong secondary flow induced by the vortex and 
depending on its intensity. The intensity of the tip vortex 
results from the hydrodynamic loading of the hydrofoil and 
on the complicated processes of vorticity concentration and 
dissipation. The results of own LDV measurements of the 
velocity field in the vicinity of tip vortices behind hydrofoils 
were used as the reference for the computations. As the result 
the computational grid structures and turbulence models best 
suited for prediction of flows dominated by the concentrated 
vorticity were selected. In the next stage the research was 
focused on prediction of the detailed geometry of the cavitating 
tip vortices. In this stage two commercial codes were used in 
parallel, namely Ansys/Fluent and Ansys/CFX. Analogically 
to the previous stage, different grid structures and different 
turbulence models were tested. The PIV measurements of 
the velocity field around the cavitating tip vortex kernel were 


used as a support of the numerical calculations. The final 
recommendations included the optimum grid structure for 
the specific task of predicting tip vortex cavitation and the 
turbulence model best suited for this task — k-e RNG. These 
recommendations constitute in fact the method for numerical 
prediction of tip vortex cavitation. 

This method is now tested in confrontation with the 
experimental observations of the tip vortex cavitation behind 
the hydrofoil. The photograph of the hydrofoil model installed 
in the measuring section of the cavitation tunnel is shown 
in Fig. 1, together with the corresponding sketch of the 
computational domain used in numerical calculations. 

Altogether nine different flow conditions were tested. They 
resulted from combination of three mean flow velocity values 
of 5.9, 5.2 and 4.3 [m/s] and three angles of attack of 4, 8 and 
12 degrees. For the lowest angle of attack only at the highest 
flow velocity some cavitation phenomena were observed, 
consequently only seven conditions are presented in detail in 
the following parts of the article. 


NUMERICAL MODEL DESCRIPTION 
— ANSYS/ FLUENT 


Calculations by Ansys/Fluent were performed using the 
unstructured grid constructed of hexahedral elements, generated 
using Hexpress Numeca. This grid consists of ~4.7 million 
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Fig. 1. Scheme of the computational domain and photograph of the 
hydrofoil in the cavitation tunnel 


STAT 
RONG TAH 


Fig. 2. Computational grid used in Fluent in the tip region of the hydrofoil 


Fig. 3. Computational grid used in Fluent 
behind the tip region of the hydrofoil 
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cells and it is refined close to the wall, so the condition y*-1 is 
fulfilled. The details of the grid in the most important tip region 
of the hydrofoil are shown in Figs 2 and 3. 

It may be seen that the density of the grid is markedly 
increased near the leading edge of the hydrofoil and in the 
region behind the hydrofoil tip, i.e. in the areas playing 
an important role in the formation of the tip vortex. The 
computations were performed using the Ansys/Fluent V12 
and employing the MUSCL (Monotone Upstream-Centered 
Scheme for Conservation Laws) scheme for convection terms 
in the transport equations. The k-e RNG turbulence model was 
used, described by the following equations: 


ô 0 
© (pk). (pku,) = 
Pa ae, Ay P ¡) 


(1) 
0 ok 
= a +G,+G, —pe- Y, +S, 
j j 
0 0 0 Os 
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Detailed description of the turbulence model and model 
constants employed in Equations (1) and (2) may be found in [ 1, 
6]. The calculations were performed for a two-phase flow using 
the mixture model. Cavitation was modeled using the Rayleigh- 
Plesset equation in Zwart-Gerber-Belamri formulation [7]. The 
gaseous phase was treated as a compressible medium using the 
perfect gas equation. 

The boundary conditions were set as follows: 


Inlet: 

- Mass flow according to the mean velocity (4.3, 5.2 or 
5.9 [m/s]) 

- Absolute temperature 283 [K] 

- Degree of turbulence 1% 

- Ratio of turbulent viscosity to laminar viscosity 10 


Outlet: 

- Static pressure 15 [kPa] (this corresponds to total pressure 
in the cavitation tunnel measuring section approximately 
30 [kPa]) 


NUMERICAL MODEL DESCRIPTION 
— ANSYS/ CFX 


Calculations by Ansys/CFX were performed using the 
unstructured grid which consisted of approximately 9 million 
elements, including about 8.2 million tetrahedral elements 
and 0.8 million prismatic elements in the boundary layer. The 
grid was generated using CFX Mesh. The calculations were 
performed using Ansys/CFX 12, making use of HRS (High 
Resolution Scheme) for the convection terms of the transport 
equations. The details of the grid in the tip region of the 
hydrofoil are shown in Figs. 4 and 5, with increased density 
in the regions of expected cavitation. 

The calculations by Ansys/CFX were performed for the 
same flow parameters and boundary conditions as by Ansys/ 
Fluent. In the case of CFX the gaseous phase was treated as 
incompressible. 


Fig. 4. Computational used in CFX in the tip region of the hydrofoil Fig. 5. Computational grid used in CFX behind the tip region 
of the hydrofoil 
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Fig. 6. Comparison of calculations by Fluent (top), CFX (middle) and experimental photograph (bottom) 
for angle of attack 4 [deg] and mean flow velocity of 5.9 [m/s] 
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COMPARISON OF NUMERICAL AND 
EXPERIMENTAL RESULTS 


In this section the comparison of numerical and experimental 
results is shown in Figs. 6 — 12 for seven flow conditions, for 
which cavitation phenomena were observed. In each Figure the 
results from Fluent and CFX calculations are shown together 
with the photograph taken in the cavitation tunnel in the 
corresponding flow condition. 

In Fig. 6 the prediction by Fluent of small region of 
cavitation area at the tip of the hydrofoil agrees reasonably 
well with experiment, while both CFX and Fluent do not 
detect the weak and detached cavitating tip vortex visible 
in the photograph. Calculation by CFX does not detect any 
cavitation in tip region of the hydrofoil. A minor root cavitation 
predicted by both programs is not confirmed by experiment. 
The reason of such difference is probably a small gap between 
the hydrofoil and the tunnel upper wall, which influences the 
flow structure locally. This root gap is not taken into account 
in the numerical model. 


In Fig 7 the calculation by Fluent agrees quite well with the 
experiment, both the small region of cavitation on the hydrofoil 
and the presence of the tip vortex are correctly predicted. 
However, the length of the cavitating tip vortex is seriously 
underestimated. The calculation by CFX underpredicts both 
hydrofoil and tip vortex cavitation. The small region of 
root cavitation indicated by CFX is again not confirmed by 
experiment. 

In Fig. 8, both numerical simulations predict quite long 
cavitating tip vortices, but they are still shorter than the one 
observed experimentally. The diameter of the cavitating tip 
vortex kernel close to the hydrofoil tip is quite well calculated. 
The prediction of the vortex by CFX is a little further from the 
experimental one than the prediction by Fluent. On the other 
hand, the prediction of hydrofoil cavitation by CFX seems to 
be closer to the observed experimentally, than it is predicted 
by Fluent. Both programs again calculate small regions of root 
cavitation, which are not confirmed by experiment. 

In Fig. 9 both programs produce practically equivalent 
predictions, which are quite close to the experimental 
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Fig. 7. Comparison of calculations by Fluent (top), CFX (middle) and experimental photograph (bottom) 
for angle of attack 8 [deg] and mean flow velocity of 4.3 [m/s] 
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observation, both in terms of hydrofoil and tip vortex cavitation. 
Similarly as before the length of the cavitating tip vortex is 
underpredicted in calculations. The diameter of the cavitating 
vortex kernel near the hydrofoil tip is quite well calculated. The 
abrupt termination of the cavitating vortex kernel as calculated 
by Fluent looks unnatural. This unnatural vortex termination 
is due to the ending of the increased density region of the 
computational grid. In case of Fluent, in order to obtain the 
converged solution within defined convergence criteria, the 
tuning of under-relaxation factors was required. The reason for 
this is the flow separation on the suction side of the hydrofoil 
and instability caused by the obtained flow structure. The 
problem is discussed further in the paper. 

As previously, the both programs indicate root cavitation, 
which cannot be seen in the photograph. 

At the highest angle of attack (12 deg), differences 
between CFD results and experimental visualization are the 
most significant. In Fig. 10, the length of the cavitating tip 
vortex by both solvers is seriously underpredicted, while its 
thickness in the region close to the hydrofoil tip is calculated 


correctly. The abrupt terminations ofa thick cavitating vortex 
in CFX as well as in Fluent look unnatural. This is again 
caused by the ending of the increased density zone of the 
computational grid and could be easily corrected if necessary. 
The radial extent of cavitation on the hydrofoil (along the 
leading eadge) is reasonable well predicted by both programs. 
In both cases, the chordwise extent of the cavitation zone is 
underpredicted. 

At the higher velocity (5.2 m/s) and angle of attack equal to 
12 deg, the cavitation on the hydrofoil and within the tip vortex 
is highly overestimated. The radial extent of cavitation on the 
hydrofoil is well predicted by both programs, but this time 
they both overestimate the chordwise extent of cavitation. On 
the other hand, the photographed cavitating tip vortex seems 
to be visibly thinner than the prediction by both programs. 
Comparing experimental visualization done for the main flow 
velocity 4.3 m/s (Fig. 10) and 5.2 m/s (Fig. 11), one can notice 
that in the second case cavitating vortex seems to be thinner, 
but less concentrated that at the lower velocity. It arises from 
the unsteady behaviour of the flow and cavitation conditions 
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Fig. 8. Comparison of calculations by Fluent (top), CFX (middle) and experimental photograph (bottom) 
for angle of attack 8 [deg] and mean flow velocity of 5.2 [m/s] 
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high sensitivity and dependence on the local flow structure 
and water aeration. 

The results for the maximum investigated hydrofoil loading 
(angle of attack 12 [deg], flow velocity 5.9 [m/s]) are not 
presented, because of lack of Fluent solver convergence caused 
by the unsteady flow separation on the suction side. 


CALCULATIONS OF UNSTEADY TIP 
VORTEX CAVITATION 


It is a well-known fact that cavitation phenomena are 
inherently unstable. This is particularly true for tip vortex 
cavitation. Although it was not possible to register the unsteady 
phenomena in the experiments, an unsteady CFD calculations 
using Fluent and CFX have been performed in order to compare 
the steady and unsteady results. The simulations have been 
performed using the previously selected k - e RNG turbulence 
model for the case of inflow velocity 5.2 [m/s] and angle of 
attack of 12 degrees. In case of both solvers the second order 
time accurate scheme was applied. The presented below results 
are obtained for time step 0.01s. 


The flow structure based on the time averaged flow field 
is shown in Fig. 12. The region of calculated flow separation 
near the root of the hydrofoil and the streamlines starting 
near the leading edge creating the tip vortex can be observed. 
Downstream of the hydrofoil, the velocity magnitude at the 
cross-section is presented. One can notice the increase of 
the flow velocity close to the tip vortex due to the secondary 
flow and the decrease of the velocity in the vortex core, as 
a consequence of the dissipation effects. The lower velocity 
areas are also shown in the wake and downstream of the 
suction side separation near the upper wall. The existence of 
separation is the reason of the convergence problem at the 
higher hydrofoil loading in Fluent simulations. The effect is 
weaker in CFX case, because the mesh is much coarser in the 
area far from the tip vortex, so the vortex structure downstream 
of the separation dissipates quickly and does not influence on 
the unsteady behaviour in this flow area. 

Anyway, it is interesting that in spite of the lack of the 
unsteady effect of the separated flow, CFX results indicate 
higher unsteadiness in the vortex area than it was obtained by 
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Fig. 9. Comparison of calculations by Fluent (top), CEX (middle) and experimental photograph (bottom) 
for angle of attack 8 [deg] and mean flow velocity of 5.9 [m/s] 
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Fluent. The fluctuations of the velocity X and Y components 
at the sections 50 mm, 200 mm and 300 mm downstream of 
the hydrofoil are shown in Figs. 13 and 14. Both velocity 
components are located in the secondary flow plane, so they 
indicate of the fluctuations in the tip vortex area. One can notice, 
that fluctuations obtained with Fluent simulations are an order 
of magnitude lower than in CFX case. The simulation time in 
CFX was 2 seconds and the fluctuations of the tip vortex with 
5 Hz frequency are noticed. Fluent simulations were done 
for much longer time ~20 seconds and no high frequency 
fluctuations are obtained in the vortex area. The unsteady effects 
are obtained close to the hydrofoil root only. 

In Figs. 13 and 14, the black line marks the area of the 
time averaged cavitating zone, i.e. the existence of the gaseous 
phase. In case of Fluent, the flow parameters distribution at the 
section (50 mm) close to the hydrofoil show the highly three 
dimensional shape of the cavitating vortex, far from the regular 
cylindrical shape. Further downstream from the hydrofoil it 
becomes “more” cylindrical. In case of CFX, such an effect 
is much weaker. 


The next difference is the maximum fluctuations location 
within the cavitating zone. In the Fluent case, the highest 
values are spread along the surface dividing the two phases. 
CFX results indicate the maximum fluctuations at the vortex 
centre. The experimental assessment of the numerical results 
accuracy can be done only by the high resolution flow velocity 
measurements within the range of cavitating vortex. Such 
measurements were not possible in the research project. 

The gaseous phase fluctuations are shown in Fig. 15. As 
regards the fluctuations distribution, the results are similar. The 
maximum values are close to the time averaged inter phase 
border surface in case of both solvers, but CFX values are 7 
times higher than obtained in Fluent. 

The time averaged secondary flow velocity is shown 
in Fig. 16. The velocity in plane normal to the main flow 
indicates the vortex intensity at the consecutive sections 
downstream of the hydrofoil. It shown that in the gaseous 
phase zone (or in fact two-phase zone) the velocity is rising 
from the vortex centre to the location where the gaseous 
phase disappears, then the velocity decreases in the water 
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Fig. 10. Comparison of calculations by Fluent (top), CFX (middle) and experimental photograph (bottom) 
for angle of attack 12 [deg] and mean flow velocity of 4.3 [m/s] 
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Fig. 11. Comparison of calculations by Fluent (top), CFX (middle) and experimental photograph (bottom) 
for angle of attack 12 [deg] and flow velocity of 5.2 [m/s] 


area. The maximum values obtained in Fluent are higher 
than in CFX. It means that the vortex intensity predicted 
by Fluent is higher, what is reflected by the cavitating zone 
extension. In Fig. 17, the iso-surface of gaseous phase for 
time averaged results is shown. It can be compared with 
Fig. 11, where the steady results are presented. In case 
of CFX, the steady and unsteady results are similar. The 
unsteady simulations results indicate a bit lower range of 
cavitation on the hydrofoil and longer cavitating zone in the 
tip vortex. In case of Fluent results, the difference is much 
higher. In spite of low fluctuations in the tip vortex area, one 
can see much smaller extension of the cavitating zone on the 
hydrofoil. It is now much more closer to the experimental 
visualization than it was obtained by steady simulations. 
The higher secondary flow intensity influences the longer 
tip vortex and the length of cavitation zone. One has to 
emphasize that three dimensional structure of the vortex and 
its skewness is also much higher than in steady simulations 
Fig. 12. Time averaged flow velocity and streamlines (Fluent) and both (steady and unsteady) CFX results. 


10 POLISH MARITIME RESEARCH, No 3/2012 


RMS X Velocity RMS Vel X 


a 0.020 
—— 0.018 
0.016 
Er 0.014 
0.012 girfa, 


“spt Coes 
E aN a 
0.006 
0.004 
0.002 
0.000 


0.08 0.1 01% ‘i 0.16 0.18 0.08 0.1 0.12 0.14 0.16 0.18 
m 


RMS X Velocity RMS Vel X 


0.020 
0.018 
0.016 
0.014 
0.012 
0.010 
0.008 
0.006 
0.004 
0.002 
0.000 


SE mee 


e SE Ene 


PT i 
e 


0.08 01 O12 014 016 018 0.08 01 O12 014 016 0.18 
X [m] X [m] 


RMS Vel X 


0.2 
0.2 
0.2 
0.1 
0.1 
0.1 
0.1 


RMS X Velocity 


0.020 
0.018 
0.016 
0.014 
0.012 
0.010 
0.008 


0.1 
0.0 
0.0 
0.0 


0.006 
0.004 
0.002 
0.000 


0.08 01 0142 014 016 018 0.08 01 O12 014 016 0.18 
X [m] X [m] 


Fig. 13. The time-dependent fluctuations of the X velocity component (Fluent — left, CFX -right) at the cross-sections (from top) 
50 mm, 200 mm and 300 mm behind the hydrofoil (black line indicates the time averaged cavitating zone) 
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Fig. 14. The time-dependent fluctuations of the Y velocity component (Fluent — left, CFX -right) at the cross-sections (from top) 
50 mm, 200 mm and 300 mm behind the hydrofoil (black line indicates the time averaged cavitating zone) 
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Fig. 15. The time-dependent fluctuations of the gaseous phase (Fluent — left, CFX -right) at the cross-sections (from top) 
50 mm, 200 mm and 300 mm behind the hydrofoil (black line indicates the time averaged cavitating zone) 
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Fig. 16. The time-averaged velocity of flow in the X-Y plane as calculated by Fluent (left) and CFX (right) in the cross-sections (from top) 
50 mm, 200 mm and 300 mm behind the hydrofoil 


14 POLISH MARITIME RESEARCH, No 3/2012 


-=— Flow direction [m] 
135 13 12 12 115 11 105 1 095 09 085 08 075 07 0.65 0.5 0.55 Lp. 


~#— Flow direction [m] 
12 115 11 105 1 095 09 085 08 075 07 065 06 055 05 
a Rs AA AAA 


Fig. 17. Gaseous phase iso-surface (time-averaged) - Fluent (top), CFX (bottom) - angle of attack 12 [deg] and flow velocity of 5.2 [m/s] 


CONCLUSIONS 


- The numerical investigations of the cavitating tip vortex 
generated on the hydrofoil are presented. It is demonstrated 
that both solvers CFX and Fluent predict the tip vortex 
and the cavitating zone qualitatively well. In both cases, 
the trend of dependence of cavitation extent and intensity 
on the inflow conditions was reflected properly. At high 
inflow velocity and high angle of attack the separation 
on the suction side of the hydrofoil was obtained, which 
is the reason of the instability and the lack of converged 
solution in Fluent case. The reason of the higher sensitivity 
to the flow structure dependence in Fluent case is the mesh 
resolution in the wake along the whole span. In CFX case, 
the mesh is refined close to the wall and in the tip vortex 
area, only. 

- The simulations for higher loading (high inflow velocity and 
angle of attack) indicate necessity of the unsteady effects 
to be taken into account if the highly three dimensional 
effects are detected in steady simulations. It leads to better 
agreement with the experimental visualization. 

- Since cavitation is highly sensitive and dependent on the 
local flow structure and water aeration, the quantitative 
comparison requires detailed information of the fluid quality 
in order to set an adequate boundary conditions for two 
phase flow. 

- It is difficult to assess the flow velocity distribution in the 
vortex as well as velocity fluctuations and gas volume 
fraction, until the high resolution in space and time flow 
field measurements are not performed. 
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ABSTRACT 


In this paper, we consider an optimization of the hull shape in order to minimize the total resistance of a ship. 
The total resistance is assumed to be the sum of the wave resistance computed on the basis of the thin-ship 
theory and the frictional resistance. Smoothness of hull lines is proved with mathematical procedure, in which 
differentials of the hull lines functions are analyzed. The wave-making resistance optimization, involving 
a genetic algorithm, uses Michell integral to calculate wave resistance. A certain hull form is generated 
by the method using cross section information of a modified DTMB model ship 5415 and a comparative 
experiment is carried out. Experimental and calculation result show that the method is of good adaptability 
for designing certain types of ships with excellent resistance performance. 


Key words: hull lines design; wave-making resistance; optimization 


INTRODUCTION 


Resistance performance has a significant effect on the 
operating cost of civil ships and the survivability of military 
ship. In the optimization design of the hull, it is the total 
resistance that people care about most [4, 19], which is 
generally the sum of the viscous resistance and the wave-making 
resistance. In a sense, once the principal dimensions of a hull are 
determined, there is no significant wetted surface change and the 
optimization design of hull form is to obtain a hull form with the 
minimum wave-making resistance [1]. Further more, for high- 
speed ship, it is reasonable to take the wave-making resistance 
as a major objective in the optimization design because of its 
high proportion in the total resistance [22]. 

Wave-making resistance can be greatly reduced by 
employing excellent hull form or optimization of some 
existing hull offsets. Shape optimization is a growing field 
of interest in many areas of research, such as marine design 
and manufacturing. Wilson, Hendrix, and Gorski [20] 
develop a computational tool set and process framework help 
designers to decide the hull shape, which had great effect on 
its hydrodynamic characteristics. 

In order to obtain a hull with the minimum wave-making 
resistance, designers often try to make improvement of parent 
hull. This requires lots of ship design experience, without 
mentioning the uncertain improvement [21]. Chen and Guang 
[2] optimized the shape of the after hull based on the desired 
wake distribution, using B-spline surface method to generate 
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the surface geometry of ship. Peri, Rossetti, and Campana [16] 
carried out numerical shape optimization of a tanker ship with 
the aid of CFD techniques and experimentally verified. The hull 
form optimization often concerns one of the most important 
applications of wave-making resistance theories. Grigoropoulos 
and Chalkias [7] developed a formal methodology for the hull 
form optimization, using parametric hull form modeling to 
generate the variant hull forms, in which Rankine-source panel 
method and strip theories were both involved. 

Development of a three-dimensional hull fairing form is 
one of the main requirements in the design of a marine vehicle. 
The final hull form must satisfy both the desired shape and 
performance characteristics, such as resistance performance 
[17]. Ghassemi and Ghiasi [6] developed a numerical program 
to determine the total resistance of planing crafts, and four 
different hull forms of Series 62 model 4666 planning craft were 
presented as calculation examples. A parametric approach to 
design of hull form was studied by Zhang, Zhu, and Leng [23], 
which provided the means for quick generation and variation 
of hull form for the hydrodynamic optimization of hull form. 
Pérez, Clemente, Suárez, and González [14] used explicit 
spline curves to make a wire model of the ship stations. An 
inverse design algorithm in determining the optimal shape of 
the bulbous bow was developed by Chen, Huang, and Fang [3], 
with the Levenberg-Marquardt method and B-spline surface 
control technique utilized. 

Genetic algorithm has been used widely in hull design 
since it appeared several yeast ago. A Multidisciplinary design 


optimization method was used to optimize the DTMB model 
ship 5415 by Peri and Campana [15], while Gammon [5] 
conducted optimization to fishing vessels applying a Multi- 
Objective Genetic Algorithm. Li [10] proposed a hybrid 
approach for multi-objective optimization of ship’s principal 
parameters in conceptual design, employing a multiple 
objective genetic algorithm. Lu, Lin, and Ji [11] calculated 
the free trim hydrostatic ship characteristics applying the 
genetic algorithm, and some necessary improvement was 
made in practice to speed up the evolution. Kim and Yang [9] 
utilized A multi-objective genetic algorithm to develop a hull 
surface modification technique for the CFD-based hull form 
optimization. 

Mathematical ship enjoys popularity among ship designers 
due to its good adaptability and excellent processing 
performance [12].The paper is going to develop a method for 
mathematical hull lines design, including the optimization 
calculation of the ship resistance. A kind of quadratic curve is 
used to generate a certain hull of the minimum wetted surface, 
and smoothness of the designed lines is discussed. A genetic 
algorithm is used to modify cross section of the hull and 
Michell integral is applied to solve the wave-making resistance 
in the optimization procedure. Numerical calculations 
and comparative experiments are conducted to assess the 
availability of the method. 


HULL LINES DESIGN 


Cross section constructed by quadratic curve 


Quadratic curve is applied to construct cross section of the 
hull; the main process is discussed in the following parts. There 
are some conditions must be meet first: 

Condition 1: The ship has only one symmetric plane called 
the centerplane. The shape of underbody of the ship is shown 
in Fig. 1. 


Fig. 1. Underbody of the hull 


Condition 2: Half of the cross section presents with certain 
types of shape as Fig. 2 shows. 


Fig. 2. Certain types of cross section constructed with arc curve 


Condition 3: Both the middle buttock line and designed 
waterline are fairing. 


A Cartesian coordinate system, as shown in Fig. 1, is defined 
with the x-axis coincident with heading direction of the ship, 
and positive distance measured upstream. Let L be hull length, 
so the value of X ranges in [0, L]. Middle line L, is indicated 
with function b(x), designed waterline L, with function a(x), 
and half of cross section area with function S(x). Within the 
scope of [0, L], a(x), b(x) and S(x) are all bounded and smooth 
functions, and their first derivatives are continuous. As the 
hull has a symmetry plane called the middle sheer plane, only 
hull lines and section areas of the starboard side (y > 0) are 
discussed. The paper is to discuss the underbody of hull, and 
lines above the designed waterline will not be involved. 

The shape of the underbody, which can be expressed by 
a bunch of body lines, is indicated with a mathematical equation 
as follows: 


z? 2 


E a 
ba} a(x)? 

Where c = c(x) is an undetermined parameter used to modify 
the cross section. A new parameter f is introduced: 


B(x) = 2S(x)/a(x)b(x) (2) 
Assuming that cross section area of the hull reaches the 


maximum when x = L/2, that is to say, there exists x, € [0, L/2] 
making B(x,) = 1. 


+czy=1 (1) 


Mathematical analysis of the approximate 
minimum wetted surface area 


To gain a hull of good resistance performance, mathematical 
process is taken to analysis character the quadratic curve that 
constructs the cross section. As we know, arc length is the 
shortest to that of any curve with the same area. A comparison 
of the length between the curves mentioned above is taken in 
this section. 

As it is shown in Fig. 2, the circular arc is defined with 
radius r and central angle 20. We have: 


rsin0 = 5 Ja} +b(x)? 


(3) 
1 1 
0 => [a(x +b(xY1cos0 + S()-5200b(%) 
Let lire be length of the circular arc, then: 
Lae a 2r0 (4) 


The quadratic curve can be expressed as: 


y(z) =—Fa(x)ez+a(x) i- e E y gp f (5) 


Let Lonc be length of the quadratic curve, then: 


b(x) 
o = [ 1 + 


Area of cross sections: 


S(x)= PE aK cz+a(x) ii- =x a4 ¿aro o Ja 3 


The discussion of Eq. (1) for x < Xo, X = Xp and x > X, is 
detailed in appendix A and Length of the two kinds of curve 
mentioned above are calculated. Solutions of the Eq. (1) are 
discussed as well and the conclusion becomes: 


conic 


[y'(z)] dz (6) 
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- when 0 < x < Xp, namely B(x) = 2S(x)/a(x)b(x) < 1, and 
c > 2/a(x)b(x) Eq. (1) has at least one nontrivial solution; 

- when x = Xo, namely: B(x) = 2S(x)/a(x)b(x) = 1, and 
c = 2/a(x)b(x) Eq. (1) has at least one nontrivial solution; 

- when x < x, < L/2, namely: B(x) = 2S(x)/a(x)b(x) > 1, and 
c < 2/a(x)b(x) Eq. (1) has at least one nontrivial solution. 


With a set of cross sections, of which breadths, draughts 
and areas are limited in a certain range, one can calculate the 
length of the designed quadratic curve and the circular arc. 
Most of the relatively error is less than 0.05% by comparing 
the designed quadratic curves with circular arcs, as illustrated 


in Fig. 3, and the maximum error is only 0.11%. 
0.15 
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Fig. 3. Length error of the designed curve to circular arc at every station 


It is not difficult to get a conclusion that hulls consisted of 
the designed cross section curves can be approximately seen 
as hulls of the minimum wetted area. 


Simple mathematical discussion of rationality of 
the designed hull lines 


In the following part, whether the designed lines are 

accordant with the practical situation is discussed, that is: 

- whether the first derivatives with respect to x of any 
waterline and buttock line are smooth. 

- whether changes of signs of the second derivative with 
respect to x of any waterline and buttock line are in 
accordance with general law of hull lines. 


A theorem for implicit differentiation is introduced firstly. 
Suppose that F(y, z) meets the following conditions: 
- in the domains |z — z| < 6 and |y — y,| < A, both F, and F, 
are continuous. 
- FY, Zo) = 0 
- Eo Zo) #0 


These lead to some conclusions: 

- ina certain neighborhood of point (yy, Zo), F(y, z) =0 uniquely 
determines a function y = f(z); and y, = f(Z,). In other words, 
y = f(z) is defined within a certain neighborhood O(z,, n) 
of Zọ meeting F(z, f(z)) = 0 and y, = f(Z); 

- y = f(z) is continuous in the domain O(z,, n); 


The first derivative of y = f(z) is continuous in the domain 
O(Z,, n), and at any point where E, # 0, y? =— F_(y, z)/F,(y, z). 
On the basis of the theorem mentioned above, the implicit 
function that c = c(x) meets is analyzed as follows: 
- when 0 < x < Xx, namely B(x) < 1, Eq. (A5) can be 
expressed as: 


F(x,c) =exp gaca + 
4 (8) 
-PORD 1 _Laggbye=0 
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- when x < x, < L/2, namely B(x) > 1, Eq. (A13) can be 
expressed as: 


Fe, nl $ ERE), 
_ 200%” _ 
¡28032 0 


For arbitrary x € [0, Xo) or x € (Xp, L/2), it is easy to validate 
that F(x, 2/a(x)b(x)) = 0, F(x, 2/a(x)b(x)) + 0, where F, and F, 
are continuous in a small enough neighborhood of x. 

In that case, it can be concluded that c = c(x) is continuously 
differentiable function in domain [0, x,) and (x,, L/2]. Since 
both equations have the solution c(x,) = 2/a(x,)b(%,), e = c(x) 
is continuous in the domain [0, L/2]. A complete process is 
detailed in appendix B. 

The first and second derivative of y can be written as: 


vn ating east Eje je 
y'(x)=-a'ah Pesa! fl E ¿Yo h 


(9) 


ah? E Hago +hatcc) (10) 
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It is clear that the waterlines of the hull are feasible from 
the discussion above. Since y and z are interchangeable in 
Eq. (1), a similar conclusion about the buttock lines can be 
got as well. 


Wave-making resistance optimization and 
sample designs 


Wave-making resistance calculation 


Since [13] published one of his scientific papers in 1898, 
the famous Michell integral has made great achievements in 
ship wave-making resistance and hull optimization. The basic 
formulas of Michell integral are written as follows [18]: 


Ape? * 
R,= or o +J’) 


TU mi (12) 
I= fín. ze" cos(Agx/U?)dxdz (13) 
J= f {n, (x, Ze el” sin(agx/U?)dxdz 04 


Where n(x, z) defines the geometry of the hull, g is 
acceleration of gravity and U indicates the moving velocity 
of hull. As have been defined above, x-axis is coincident with 
heading direction of the ship, and z-axis is coincident with 
draught direction the ship. 

In the research, a numerical program is conducted applying 
Michell integral to calculate wave-making resistance of hulls. 
Before the calculation, the Wigley N43 hull is used to check 
the accuracy of the numerical procedure. Line function of the 
hull can be written as [8]: 


las fo fred), 
(2) af) 


Wave-making resistance coefficient C, is defined as: 


w 
E o DV? 

where R, is wave-making resistance of the hull in N, D is the 
hull displacement in kg, and V is the heading velocity of the 
ship in m/s. Curves of C,, versus Froude Number are presented 
in Fig. 4 to compare the result calculated in the program with 
that Wigley got from the calculations and experiments, with the 
ordinate representing the coefficient of wave-making resistance, 
and the abscissa the Froude number. 
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Fig. 4. Wave-making resistance comparison of present calculation with that 
of Wigley ss and the experimental result for Wigley N43 model 


As it is illustrated in Fig. 4, though the agreement of the 
calculation at low Froude Number is not so perfect, the result 
show in this paper shows good agreement with that Wigley got 
at high Froude Number (Fr > 0.32). 


Genetic algorithm 


A genetic algorithm optimization model is applied and the 
numerical procedure is shown as follows: 

The hull is divided into 21 stations in the direction of length. 
Since the section area of the No.0 station and No.20 station are 
zero, section areas of No.1 station to No.19 station are defined 
as designed variables: 


X=(11),1=1,Z 319 (17) 


Michell integration is used as the target function, that is, 
the adaptive function of genetic algorithm, and the volume 
of displacement of the ship as the constraint condition. 
5 percentages more of section area of parent hull is taken as 
the upper limit of the optimum region (X,,) while 5 percentages 
less of that as the lower limit of the optimum region (X, ).Then, 
the optimum model can be expressed as: 


Target funktion Min R,=R,(X, U) 
Subject to gi X> A, (18) 
g: Xi < Xy 


gi: (A — A/A < 5%o 


Where Rw is wave-making resistance, U is the moving 
velocity of hull and A indicates the volume displacement of the 
hull. Since the volume of displacement constraint condition is 
nonlinear, we use the punishment coefficient method to simplify 
the computational process, and the punishment function can 
be written as 


AA, pl 
Rw(X, U, A)=Rw+M x (19) 
Where M and N are coefficient that need to be determined. 
Since it is difficult for genetic algorithm to deal with both 
M and N, we can make N constant and M a kind of descending 
series. 


Design sample and analysis 


Lines designing example of a hull is shown in the following, 
with the modified DTMB model ship 5415 used as the parent 
hull, of which the bulbous bow is removed. The main principles 
are listed in Table 1. 


Tab. 1. Main parameters of the parent ship 


Parameter 


Displacement [t] 


12129 


Length overall [m] 


172.5 


Beam [m] 


23.09 


Depth [m] 


12.34 


Volume of displacement [m°] 


11785.6 


Amid ship section coefficient C,, 


0.817 


Length on the waterline [m] 


160.0 


Waterline breadth [m] 


21.46 


Draught [m] 


6.927 


Prismatic coefficient C, 


Block coefficient C, 0.496 
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0.607 
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Since the hull lines of new ship generated is similar to the 
parent hull as it is shown in Fig. 5, resistance test of 1/30 scale 
models of the two ships are performed in the towing tank of 
Huazhong University of Science & Technology respectively for 
further comparison and analysis. The experimental temperature 
is about 18.0°C, and the models are kept afloat on an even 
keel during the test. Some useful parameters are shown in 
Table 2. 


Parent Ship 


Fig. 5. Lines ian the parent nn and the designed one 


Tab. 2. Some useful information for the experiment 


Parent ship 
model 


New ship 


Parameter model 


Wetted surface area [m?] 


Displacement [kg] 


Average draught [m] 


The resistance of the models tested in the experiment is 
converted to that of the full-scale ships, and the total resistance 
coefficient and residual resistance coefficient of the hulls are 
shown in Fig. 6. Experimental result indicates that the resistance 
of the new ship is similar to that of the parent hull. With the 
velocity of full-scale ship ranging from 22 to 36 knots, the total 
resistance of designed hull is 1% averagely more than that of the 
parent ship. Taking the average accuracy of model experiment 
into account, a conclusion can be drawn that resistance of the 
designed ship has no obvious change compared with that of 
the parent hull (Fig. 6). Obviously, the wetted surface area 
of the new ship declines slightly. These is probably because: 
firstly, the wetted surface area coefficient of the parent ship 
model is too small to bear too much change; secondly, on the 
condition that the displacement and the principal dimensions 
are almost the same, when the character of the body lines 
changed a little, it leads to no big difference in decreasing the 
wetted surface area. 
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Fig. 6. Comparison of resistance coefficient of the new hull 
with that of the parent one 
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Optimization of designed ship of full-scale is conducted, 
using Michell integral and the genetic algorithm mentioned 
above. Some different velocities such as 26kn (Fr = 0.338), 
28kn (Fr = 0.364), 30kn (Fr = 0.390), 32kn (Fr = 0.416) are 
taken into account. As have been stated, the wave-making 
resistance is taken as the target function. New section area 
curves of different resistance performance are obtained. The 
cross section area curves of different velocities are shown in 
Fig. 7, compared with that of parent ship. 
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Fig. 7. Transverse-section area curves of different velocities optimization 
compared with that of parent ship 


Comparison wave-making resistance coefficient is shown 
in Fig. 8. The calculation result of 26kn shows good resistance 
performance at low velocity, while the wave-making resistance 
at high speed has no obvious improvement. Calculation result 
of 28kn and 30kn show just a little change at low speed in 
wave-making resistance, but the improvement at high speed 
is noticeable. Based on an overall consideration of different 
velocities, optimization result of 28kn is chosen for hull lines 
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Fig. 8. Coefficient of wave-making resistance of different velocities 
optimization compared with that of parent ship 


Wave-making resistance coefficients of parent hull and 
the optimized one at 28 knots, compared with experimental 
result, are shown in Fig. 9. Lines of the designed hull and the 
optimized one are shown in Fig. 10. Since cross section areas 
of the hull vary in a small range, lines of hull dose not have 


much change. That is to say, with just a little change of hull 
lines, the hull can get some noticeable improvement in wave 
making resistance performance. 
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Fig. 9. Comparison of calculation result with experimental result 
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Fig. 10. Lines of the designed hull and the optimized one 


In general, the wave-making resistance performance of 
the new hull is close to that of the parent hull. At low Froude 
number, the wave-making resistance coefficient of optimized 
hull is a little smaller. However, when the Froude number is 
higher, the coefficient of optimized hull is apparently smaller 
than that of parent one. That is to say the wave-making 
resistance performance of optimized hull has been obviously 
improved though the result coincides well with experimental 
result at low Froude number but not very well at high Froude 
number. 


CONCLUSION 


A hull lines design method, which involved wetted surface 
analysis and wave-making optimization, is applied with the 
intent of designing hulls of good resistance performance 
using a given set of constraints. A certain kind of quadratic 
curve is used to generate hull lines, of which the feasibility is 
proved in a brief theoretical analysis; approximate minimum 
wetted surface of the hull is calculated as well. Hull lines are 
proved smooth through mathematical procedure, in which 
both the first and second differentials of the line functions 
are derived. A program is compiled to calculate wave-making 
resistance of hulls using Michell integral and the Wigley hull 
model is taken as an example of validation. ADTMB model 
ship 5415 is used as an example of hull form optimization 
designing, applying a genetic algorithm. The particular genetic 
algorithm developed during this study uses the wave-making 
resistance as the objective function and allows automatic 
modification of cross section curve. Experiment result of new 
hull and the parent hull indicate that the hull optimization 
designing method is practical. With respect to the calculation 
result, the wave-making optimized hull is of good resistance 
performance. 


Appendix A 


Discussion of Eq. (1) for x < Xo, X = x, and x > x, is shown as follows: 
1) When x < x,, namely B(x) = 2S(x)/a(x)b(x) < 1: 


- when c < 2/a(x)b(x) namely 1/b(x)? — 


$i E.. E 


Since there exists a unique zero solution of equation z = sinfz, then 1/b(x)? — 


216) ares 


sol . | aE SEG 


a(x)?c2/4 > 0: 
1—a(x)’b(x)’c? 
4 


contradicts the previous assumptions. Thus, Eq. (1) has no nontrivial solution. 


- when c = 2/a(x)b(x) namely 1/b(x)? — 


then B(x) = 1 and it contradicts the previous assumptions. Eq. (1) has no nontrivial solution. 


- when c > 2/a(x)b(x), namely 1/b(x)? — a(x)202/4 < 0: 


sy- 


-arcsin (Al) 

then: 
(A2) 

a(x)?c?/4 = 0 i.e. c = 2/a(x)b(x), which 
a(x)2c?/4 = 0: 
S(x) = a(x)b(x)/2 (A3) 
labe e , 1 

nl e ES 1+ > a(x) b(x)c (A4) 


[a bay E 
2 BE i. 
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then: 


exp(B POBA y EA +a(x)b(x)c (AS) 


By using Taylor expansion, an approximate solution of c can be got: 
3-3B —|—-15B (x)* +248 &)*-18B 69” +9 ¡ 
+1 
2(B)x* 
a(x)b(x) 


Since the equation e® = z+-/z? +1 has nontrivial solutions when f < 1, the conclusion comes out that Eq. (1) has solutions. 
Loonie can be solute by considering Eq. (6) and (A6). 


conic 


2 


c(x)= (A6) 


2) When x = Xp, namely B(x) = 2S(x)/a(x)b(x) = 1: 
- when c < 2/a(x)b(x), namely 1/b(x)? — a(x)2c2/4 > 0 


ak) cates 1—a(x)*b(x)?c? 


eGo a(xy 0/4 4 


then: 


E h- ats) eye 7 pera (A8) 


Since there exists a unique zero solution of equation z = sin z, then 1/b(x)? — a(x)?c2/4 = 0 Le. c = 2/a(x)b(x), which 
contradicts the previous assumptions. Eq. (1) has no nontrivial solution. 


(A7) 


- when c = 2/a(x)b(x), namely 1/b(x)2 — a(x)202/4 = 0. It is easy to verify that Eq. (1) has at least one nontrivial solution. 
Substituting the equation above into Eq. (1): 


EIA q yZ= (A9) 
xy ACK)? ACK) b(Ko) 
namely: 
2 
ae y j] =] (A10) 
b(%) a(xo) 


The designed curve and circular arc are both turned to lines of the same expression z/b(x,) + y/a(x,) = 1, and the curve length 


koloa =V +b. 


conic — “circle 


- when c > 2/a(x)b(x), namely 1/b(x)? — a(x)202/4 < 0: 


__ bD | a) bec? 1 (A11) 
S(x) s O = n 4 1+ T 
cortz 


then: 


exp (B Cako = [RO B e y Faba (A12) 


Since there exists a unique zero solution of equation e7 = z+4 z? +1, then 1/b(x)? — a(x)20Y/4 = 0 i.e. c = 2/a(x)b(x), which 
contradicts the previous assumptions. Thus, Eq. (1) has no nontrivial solution. 


3) When x, < x, namely B(x) = 2S(x)/a(x)b(x) > 1 
- when c < 2/a(x)b(x), namely 1/b(x)? — a(x)202/4 > 0: 


Aa _ l- Lay bc” (A13) 


“5 [Ube alae? Poy" a(x)’ c? 


then: 


sin b qli- abe E =,/1- ater bare a (A14) 


By using Taylor expansion, an approximate solution of c can be got: 
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; 20969 —,/120B(x) — 20B(xy 
e B(x) 


a(x)b(x) 


Since the equation z = sinz has nontrivial solutions when B > 1, Eq. (1) has at least a nontrivial solution. Curve length 1 


be calculated by solving Eq. (6) and (A15). 
- when c = 2/a(x)b(x), namely 1/b(x)? — a(x)2c2/4 = 0: 


S(x) = a(x)b(x)/2 


That is B(x) = 1, which contradicts the previous assumption B(x) > 1. Thus, Eq. (1) has no nontrivial solution. 


- when c > 2/a(x)b(x), namely 1/b(x)? — a(x)202/4 < 0: 


____a(x)b() EOS ll 
S(x) MECO =n i 1+ sabe 
4 


then: 


exp (p, 800 BOY y) =, OLEOLE 


Since the equation e = z+-~/z? +1 has no solutions when ß > 1, Eq. (1) has no solutions yet. 
Appendix B 


Discussion of the continuity of waterlines function is shown as follows: 
The first derivative of c with respect to z is calculated as follows: 
1) When x e [0, xo), c = c(x) is defined by Eq. (8). Let: 


O®=, A aK bR e)? -1 


Eq. (8) can be written as: 


Bx) O®) =In(O(x) +06) +1) 


Differentiate both sides of the equation with respect to x, and then we have: 
B'EJOC)YOG) + 1 


1- B(x) JOY +1 


Let limO(x)O'(x) = I. A simple form of Eq. (B2) can be got by using L’Hôpital’s Rule: 


O(x)O'(x)= 


21 
I=lim O(x)O'(x) = B'E) === 
oer (x) (x) B ( o) _ B'&,) 2j 
namely: 
lim O(x)O'(x) = I = -3B'(Xo) 
then: 


aja) A + a(x)"b(x)b'(x)e(x)’ + a(x)"b(x)'e(x)e'(x) 


O(x) = 
4,11 — 4 a(x) DE) Ax) 


Solve Eq. (B1) and (B5), then: 
va FOO (x) — a'(x)a(x)b(x)'¢(x)’ — a LD rar 
aa a(x) DEJA) 


Simplify the equation, and we can get the left-hand limit of c’(x) at xo: 
' +a! + ' 
de OB'Ko) + a'b ex) + aC%p)b' pep) 
xxi a(x, )b(X) 
2) When x € (Xp, L/2], c = c(x) is defined by Eq. (9). Let: 


ow =4,/1 4 aba) 
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conic 


(A15) 


can 


(A16) 


(A17) 


(A18) 


(B1) 


(B2) 


(B3) 


(B4) 


(B5) 


(B6) 


(B7) 


(B8) 
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Eq. (9) can be written as: 


sin (B(x)O(x)) = Ox) (B9) 
Differentiate both sides of the equation with respect to z, and then we have: 
'(x)O(x) cos(B(x)O(x 
00001) = BOL c0s(B(HIO6) (B10) 


1 — BE9cos(B(x)0()) 


Let limO(x)O'(x) = I. A simple form of Eq. (B9) can be got by using L’Hopital’s Rule: 


1=lim 069009 = BOS) RG B1) 
namely: 
lim O(x)O'(x) = I = 3f'(X,) (B12) 
° then: 


_al(x)a(x)b(x)'e(x)’ + a(x)b(x)b'(X)e(x)" + a(x)"b(x)'e(x)e'(x) 


O'(x) = 
4/1 Lacy bay) 


Solve Eq. (B8) and (B13): 
40(x)O'(x) + a'(x)a(x)b(x)'e(x) + a(x)"b(x)b'(x)e(x)” 
GG o —==————————_ O SS 
a(x) b(x) o(x) 
Simplify the equation, and we can get the right-hand limit of c’(x) at x9: 
lim e = -EPE + ADA) + abee) 
aj axb) 


From Eq. (B7) and (B14), we have lim c'(x) = lim c'(x), namely the left-hand and right-hand limits are equal as x — Xp, that is 
x>X XX} 


(B13) 


(B14) 


to say c’(x) is continuous in domain [Xp, L/2]. 
When z = h, the waterline defined by Eq. (1) can be expressed as: 


ee aU Y , x € [0, L/2] (B15) 


y(x) = tao) Hi 


Since a(x), b(x), c(x) and their first derivatives are all continuous, and it is not difficult to prove that the second derivatives of 


a(x), b(x), c(x) is continuous, the first and second derivative of y can be written as: 
Iis La 1_1 xe Ih? 
y'&)=-a'ah— zx c'h+a' 1- Jr. 
yy” ef E 
a2 o) 
(+ 
2,1 


" a. Br _ hate. " 1 " 9l 1 2 
y"(x)=—a" cz—2aa'c'z—aca"z ¿40 z+a l- E a 


at (B+ jaate?+ jee +d nate? e 


[0, L/2] (B16) 


a Mae 2 AMAT 2 a Bim 
(tae ja (1 Jae) 
1 1 6b” _ 2b" 1 2 tt i " 1 ' 1 " 
E yA c’ — 2aca'e -500a -732e ¿Yoo le , x € [0, L/2] 
1-| ate? |z? 
b 4 
Both Eq. (B16) and (B17) are smooth and continuous. 
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ABSTRACT 


The aim of this paper is to provide a simple method for the definition of a rudder for tuna purse seiners. 

The model achieved by the application of this method ensures the manoeuvrability of initial turning, 

course keeping and yaw checking. For this, the results obtained by other authors, through rudder tests and 

manoeuvrability tests of fishing vessels have been integrated with the recommendations of the International 

Normative. Finally, in the paper the method is applied to the case of a tuna vessel, providing a rudder model 

which could be optimised with CFD (Computer Fluid Dynamics) in further tests. The method proposed 
can be applied to other vessels whose main dimensions are met. 


Key words: rudder optimisation; vessel manoeuvrability; tuna vessel design 


INTRODUCTION 


Limited manoeuvrability is one of the most significant weak 
points in the operation of the fishing vessels when compared to 
other vessels. This is mainly due to the activity of fishing, which 
means shipping with nets. This characteristic not only leads to 
inefficient operation and therefore to high fuel consumption, 
but also to a higher risk of accidents. 

Above all, these threats stand out in the case of tuna purse 
seiners because of their particular manoeuvrability needs 
throughout their entire voyage; during navigation (the rapid 
pursuit of fish shoals), during manoeuvres (the releasing of nets 
at high speed) and whilst tacking (the stability of the vessel 
during collection of catch). 

In order to minimise these weak areas, the current 
construction tendency in these types of vessel is based on 
wider dimensions, optimised hydrodynamic behaviour, usage 
of controllable propellers and high-performance rudders. The 
objective of these improvements is to reduce running costs and 
to increase safety during operation. 

The main purpose of this paper is to contribute to the 
definition of high-performance rudders for tuna fishing 
vessels. For that, a simple definition method will be shown of 
an optimised rudder model for this type of boat, which serves 
as a starter point for further trials with CFDs. 

The method not only takes into account the manoeuvrability 
recommendations provided by the International Rules (IMO, 
SOLAS, and DNV), but also the results achieved in previous 
experimental manoeuvrability test studies with fishing vessels 
and rudders. 
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As aresult of the application of this method, the first step 
obtains the minimum forces necessary of the rudder to meet 
certain existing manoeuvrability requirements. Then, the data 
from previous studies will be taken into account in order to 
relate these forces with the rudder geometric characteristics 
and its relative position in the vessel. 

The rudder model obtained ensures meeting certain 
manoeuvrability demands and the Classification Societies’ (CS) 
requirements. However, other requirements are not ensured by 
the model, for this reason some requirements must be checked 
through testing with CFDs. 

This paper shows the base application method for a tuna 
fishing vessel, thus obtaining the definition of a rudder 
model. 

The method presented can easily be repeated to define 
rudders for tuna fishing vessels, guaranteeing their compliance 
to certain operational characteristics. 


METHODS 


The manoeuvrability requirements considered have 
been determined by SOLAS Part C [7], IMO MSC/ 
Circ.1053 and MSC.137 (76) [5, 6] and the Regulations of 
the Classification Society Det Norske Veritas [4] (Part 3, 
Chapter 3, Section 2). In addition to the previous compulsory 
regulations and recommendations, the results achieved 
during manoeuvrability tests in fishing vessels carried out 
in the “El Pardo” model basin [2, 3] have been considered, 
together with the conclusions extracted from other rudder 
tests [8, 9]. 


As a consequence of the integrated recommendations 
provided by these sources, a series of useful expressions will 
be obtained. These expressions will allow the definition of the 
rudder geometry, and its relative position in the vessel ensures 
the realization of certain manoeuvrability requirements such 
as: 

Turning ability. Itis a critical issue for tuna purse seiners due 
to their operational activity. Adhering to the Gertler acceptability 
criteria, the following expression can be used [2]: 


D < (-5C, + 7.2) * Ly, (1) 


The turning circle diameter value can be calculated through 
the expression [3]: 


25 
Z(1+ (Top pp S i) 0) 
D _ 0.048Lp Cy pP Pe 


Lp  sen(2a)B` (En) Lo, 


The initial capacity for manoeuvrability is regulated by 
the IMO [7]. In accordance with this, the advance (A,) must 
not exceed 4.5 times the length of the vessel, and the tactical 
diameter (D,) must be less than 5 times the turning circle 
length. 

Therefore, as a result of model trials and sea testing, the 
following expression is found: 


D, = 1.65D C, + 0.08 6) 


After analysis of turn manoeuvre tests with models [2] 
for boats with block coefficients of less than 0.6 (which is the 
case of tuna fishing vessels), the following relations can be 
concluded: 


D, =D, +0.55D (4) 
D,=0.5A, (5) 


This way one can define the advance as a function of the 
turning diameter (D). Taking into consideration the limitations 
of the D/Lpp ratio and the advance value limit in relation to 
vessel length, the maximum value turning circle diameter value 
can now be expressed, as well as the normal minimum force 
per angle unit. 

Course keeping. This can be defined as the ability to 
maintain a selected straight line course. By varying the rudder 
degree angle (spiral test) the progression can be observed in 
relation to the changes. It is desirable that the evolution is 
stable, positive and that hysteresis does not occur. The degree 
of leeway for this quality can be determined via the evaluation 
of the hysteresis loop width (a). 


46.43 
a = 18.12 - > (6) 
T 
Where: 
, 1 
T= d+) (7) 
a 
and: 
Ft 


Crt — a 
a sp-Ar-Vi 6) 
In order to minimise this value it is necessary to maximise 
(Ft/a). 
The yaw checking ability. This ability is verified by the 
zigzag manoeuvre test, in which moderate changes of course 
within time and space are measured. The initial zigzag 


manoeuvre can be evaluated by measuring the number P of 
Norrbin, which defines the angle of course turned per unit of 
rudder angle used, once a determined length is navigated [3]: 


P=K'(1-T'+T'e-T) > 0.275 (9) 


The number P is also accepted by IMO [6] to evaluate the 
initial progression capacity for manoeuvrability. 

Through the Nomoto equation, the following expression 
is obtained [3]: 


Lpp 1 


Y2 1 y2 10 
K?+K?  ® 

This defines the minimum value for the lift force on the 
rudder by angle unit (Ft/a) 

In order to ensure these forces per angle unit of the operation 
of the rudder (see Fig. 1) itis necessary to relate the said forces 
to design parameters that are controllable: the relative position 
of the rudder in the vessel, and the dimensional geometric 
characteristics of the rudder. 


a = angle of attack 


Fig. 1. Notation of forces on the rudder 


Rudder behaviour is determined to a great extent by the 
conditions in which it operates, and these are defined greatly 
by the relative position of the rudder in the vessel, the propeller 
and the hull [2, 9, 11]. 

The principle objectives aimed for in the definition of 
a high performance rudder are the following: the increase of 
lift force and the decrease of drag force. This resistance is 
greatly conditioned by the boundary layer of water which is 
created around the rudder surface. If there is turbulent water 
flow (high Re), the profile velocity increases and the pressure 
falls (Bernoulli), which is clearly unfavourable as the drag force 
increases and lift force decreases. 


Re =c.Vr/v (11) 


Another motive for avoiding turbulent water flow is the 
appearance of cavitation. This phenomenon is determined by 
the profile pressures of the rudder. Profile pressure varies with 
the rudder section, and in addition to this cavitation can be 
found in the propeller wake caused by water accelerating from 
the root to the extremes of the blades. For this we encounter 
the lowest flow pressures at the extremes of the propeller, and 
therefore the highest risk of cavitation. Sections of the rudder 
that coincide with the extremes of the propeller will be critical 
from the point of view of cavitation. 

The influence of rudder angles in cavitation is noticeable, 
with the risk of cavitation increasing the greater the rudder 
angle. Rudder angles have an important bearing on tuna purse 
seiners and the risk of cavitation should be avoided. 


o <-C, (12) 
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_Par+p-8:hg=P, _ Po— Py 


0.5pV2 ~ 0.5pV2 (13) 
P = (0.5pV2) ” 


The relation between pressure and lift coefficients can be 
expressed in the function of distance (+) from the origin to the 
profile surface by: 


c 
sj | (=) dl (15) 


Furthermore, the effect of the propeller on rudder behaviour 
is significant. The propeller creates axial thrust to the water, 
so the effect is that the water arriving at the rudder has higher 
velocity. Using Bernoulli and the Gutsche correction, we arrive 
at the following expression [9]: 


V, = V(1 — 0.8C, + 0.26): 


16 
1+0.5 2 wa a 
(1+05+ 55 (1 + =a ie 
14 
d 
V(1 — 0.8Cb + 0.26 
(ee (17) 


Nd 


In order to try to minimise the effects of turbulent flows, 
it is necessary to achieve a low V, value. For this a low x/d 
is recommendable, taking into account the minimum values 
given by CS which must be fulfilled to avoid problems of 
vibration. 

Tuna vessels however achieve high J values (0.35 to 0.94). 
For these ranges of J, experimental rudder tests [8, 9] have 
demonstrated that the value of C,/a versus the relation X/d 
increases proportionally (C,/a decreasing) to X/d = 0.4 (for 
low values of J the trend is reversed). 

Considering tuna vessel characteristics (Table 1), and the 
previous points, the initial value selected is X/d = 0.22, which 
complies with the minimum requirements of CS, and would 
avoid reaching high Re without the penalising lift force. 

The flow straightening effect also has to be taken into 
account [10], influenced greatly by the stern of the vessel, 
which causes a reduction in water flow speed when reaching 
the propeller (see Fig. 2). The usual consequence of the 
straightening effect is to increase the rudder attack angle (a,). 


Its effect can be measured by using the flow straightening factor 
(y), which relates the ship drift angle, Bg, with the flow attack 
angle to the ship axis for zero lift (aç). 


y = ay Pr 


a, =a—Y-Br 


(18) 


(19) 


As can be seen, the value of (y) should be as small as 
possible to widen a.. 

The y values obtained during trials, carried out with rudders 
that protrude vertically [10] in the superior extreme of the 
propeller diameter, have been high. This leads to the election 
of the Z/d value that indicates the relative position between the 
rudder and propeller. According to the trial results, it would be 
advisable to place the rudder so that it does not protrude beyond 
the highest blade tip of the propeller (Z/d > 1) to obtain a low 
flow straightening factor. However, it is recommended that the 
ratio Z/d < 1 is used in order to guarantee the coverage of the 
rudder span. As a consequence to the previous points, the ratio 
Z/d = 1 can also be recommended when considering bending 
movements and rudder blade height. Finally, the ratio Y/d = 0 
obtained from rudder tests [10] has achieved the highest values 
of C,/a, and the minimum values for a, and therefore also the 
minimum values of y. 

In order to determine value limits for forces on the rudder, 
it will be necessary to act on its geometry and to understand 
the consequences of this. For this, a real tuna purse seiner will 
be used an example (Table 1). 


Tab. 1. Characteristics of the tuna purse seiner “Draco” 
(Source: Infomarine, May 2006) 


X, (m) 
Main engine Power 


-2.55 
6000 kW (750 rpm) 
Controllable 152 rpm (4 blades) 


Propeller 


Propeller Diameter (m) 4.3 


The taper ratio (TR) is the relation between the lower and 
higher chords. If the value of this coefficient increases, both 
the lift force and the separation angle of water flow also rise. 


/ Inflow 
< 


Fig. 2. Angles of the flow straightening effect 
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This also leads to an increase in the bending moment, and the 
majority of the rudder is adversely affected by the propeller 
wake. Due to these reasons the ratio TR = 1 has been taken 
for the calculations. 

Another characteristic to be defined is the leading edges of 
the rudder. a square leading edge shape has been used for high 
rudder angles, necessary for the operation of tuna vessels; the 
drag force is less than for rudders with rounded edges. 

Finally, a rudder with a flap has been selected due to the 
quick turning speed and rapid evolution, which makes it ideal 
for the high manoeuvrability needed by tuna purse seiners. 
Furthermore, the rudder is going to support minimum vertical 
forces and bending moments. 


Calculation of the minimum normal force per 
rudder angle unit 


This force is calculated bearing in mind the turning circle 
abilities applicable to the acceptability criteria of Gertler, 
Eq. (1), to the vessel base (Table 1). 


D/L<-5C,+7.2=4.5m 


Therefore, D,,,. = 372.15 m is obtained. With this value, 
for a turn angle of 35° (0.61 rad), Eq. (2) must be applied. For 
this, using the vessel base, a minimum normal force value is 
obtained of: 


dF /da = 1545.31 kN/rad 


Alternatively, to verify this maximum diameter obtained 
complies with the IMO requirements regarding turning 
ability: 

A, <4.5L = 372.15 m 


Now considering expressions (3) — (5) and substituting the 
values, the result is: 


Dox = 372.15 < 601.76 m 


Therefore, for the minimum value of F /a would lead to the 
turning ability requirements of the IMO being met. 

To calculate the minimum lift per angle unit for the correct 
operation of the vessel, firstly the number P through Eqs. (9) 
and (10) is analysed. Nevertheless, in order to find the minimum 
value of lift force by turn angle, it is necessary to define the 
relations between K’ & T’ with F/a. To obtain these relations, 
expressions (16) and (17) are integrated into the following 
formulae: 

K? +K? 1 
—— (20) 


0.27+0.258(F+0.38)0390.54 pe 6+1 


pp 


T= 


G = 0.019324 (1 a a 
pp 


Fr : 
a _ (1.26 — 0.8C,)?- 


F = T 


1 
7 P Ve 


1\/ 1 25(Tpp — Tpr) 
(Gt 
B/ \Cp Lpp 
The previous expressions have been empirically obtained 
from reference [3]. 


Considering also the values of the vessel base (Table 1) the 
relations are finally obtained: 


(22) 


0.065 


r= 0.39 (23) 


0.258(7.19 : 10-64 + 0.38) 0.26 


K’ = 3.788 - 1077. 
Fe 
Qa 


| 0.258 (7.19 - 10-6 Fe + 0.38) 


(24) 
0.39 


— 0.22 


The P number can be expressed as a function of F/a. In 
order to evaluate the lower limit of this value, the functions 
which make up expression (9) will be analysed: 


P=(1-T’+T’e")K’ =f, - f, > 0.275 
Where: 
f,=1-T’+T’e!7 (25) 
and: 
f, = K’ (26) 
Function f, (Fig. 3) is not continuous on T’ = 0 and its limits 


tend infinitely towards 0 and 1. 


T 
Fig. 3. Tendency of function f, of the number P of Norrbin against T” 
If function f, (Fig. 4) is analysed, which depends on F/a, it 


can be seen that a discontinuity in f, = K’ = 0 also exists. The 
maximum of function f, is obtained by F/a = 441 kN/rad. 


400000 600000 800000 1000000 


Ft/a (N/rad) 


Fig. 4. Tendency of function f, of the number P 
of Norrbin as a function of F/a 


This value of F /a signifies that T’ < 0 and, therefore f, takes 
negative values. To ensure T’ > 0 it will be necessary to use 
F/a > 838.3 kN/rad, even though this value does not comply 
with the requisite given in Eq. (9). For this it is necessary to use 
expression F/a > 2009 kN/rad (Fig. 5), which will correspond 
to the minimum lift force required per unit angle. 
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0 1000000 2000000 3000000 


Ft/a(N/rad) 


4000000 5000000 6000000 


Fig. 5. Tendency of number P of Norrbin against Ft/a (N/rad) 


Although the increase in the F/a value is favourable for 
complying with expression (9), this increase penalizes the 
ability to keep course (6). In addition, these high values of F/a 
implicate an enlargement of the servomotor and rudder. 

The leeway for the course keeping ability is measured by 
the width of the hysteresis loop, which only has meaning for 
positive values. Despite the fact that expression (6) is zero for 
T’ = 2.56 (Fig. 4), for this value expression (9) (P = 0.128, 
fi = 0.174 y f, = 0.737) is not carried out. 


T 


Fig. 6. Width tendency of the hysteresis loop against T? 


Following the previous results, the minimum lift force per 
angle turn was assumed by: 


(F/O) min = 2744 KN/rad 


This allows the value to comply with expression (9) without 
excessive risk to the course keeping ability (6). 


Calculation of main characteristics 


The fundamental geometric characteristics that identify 
a rudder are: span (h) defined in the normal flow direction, 
chord (c) which is the measurement of the rudder blade, 
thickness (t), perpendicular to the longitudinal axis of the vessel 
(see Fig. 7). Other parameters are: the profile type (t/c), area of 
the rudder (A,), defined as the product span by the chord: 


A, =h.c (27) 
and the aspect ratio (A), which is the relation between the rudder 
span and the average chord measurement (elongation): 

X= h/c (28) 
Firstly, in order to select the rudder characteristics it is 
necessary to take into account the following points: 
+ The stern post characteristics (for a relation ratio Z/d = 1) 
limit the span of the rudder (h) to 6 m. 


¢ In accordance with CS, the minimum area of the rudder 
blade [5] is limited to A, = 8.26 m? for the vessel base 


according to: 
2 B ; 
1 + 50C2 | — 
Lpp 
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TL (29) 
A. = —PP 
j 100 


centerline stock 


Fig. 7. Notation of the main parameters in the rudder 


e [tis recommended that the elongation ratio (A) is as high as 
possible because of the increase in the lift force. However, 
on the other hand, for flap rudders the chord must be wider. 
Due to the restriction of the space on tuna vessels, the 
possible maximum span is 6 m. 

e Observing the flow speed on the rudder and the possible 
chords (16), the Re values are expected to be high. 
Therefore, as the lift coefficient increases (C,,) the drag 
decreases (C,), and the separation angle (a,) and cavitation 
risk also increase (see Table 2). 

e Finally, note that for equal values of Re, the drag force 
coefficient (Ca) increases when the t/c ratio also rises [11]. 


Fa 
o Ya 
FPA - VŽ 


Table 2 shows different calculated possibilities of rudders 
for different Re values (and therefore chords), analysing their 
respective stall angles. These have been calculated through 
expressions (17), (27) and (28) integrated into [3]: 

1.2 =) 


t 
as(degrees) = 7.11 (1 +7. -) (1 + ZS) 


8K 0.5 h 
1+ 0.048 (1 (1 + =) = 


The previous expression has been obtained from model 
tests and verified with real fishing vessels up to 6 m of rudder 
span. In addition to this, for the t/c calculation expressions (8) 
and (16) have been integrated together with the minimum lift 
force value obtained in [3]: 


“tay, = 5 (30) 


(31) 


Ft 
Ce a 


= a = 0.357). 


a ZApVE  2+255 
Cp +0.3 


8 
(mJ?) 
(14K h -0(5) 
14121910 y ) 


Finally, the expression used is: 


(32) 


À +2.55 


= =a oF 
(1-12.41 Tm) 


Tab. 2. Rudder geometric characteristics for different Re values with their 
respective stall angles 


Table 2 shows the valid solutions achieved which meet the 
maximum span rudder allowed (6 m) with relation t/c remaining 
positive. The chosen option is that highlighted in Tab. 2. This 
option is the most interesting for a flap rudder as it has the wider 
chord (c), a not very high t/c profile relation and the highest 
stall angle (critical for tuna vessels). This option also complies 
with the minimum rudder area demanded by CS. 


Stall angle correction 


The stall angles, indicated previously, have been calculated 
using empirical expressions obtained in successive tests on tuna 
vessels of fixed rudders blades with taper ratios different to 
one. Therefore, one must consider the effect of the rudder with 
a taper ratio value equal to one in the stall angle calculation. 
For this, the following expressions (suitable for rudders with 
a squared blade tip) have been used [9], where the angles shown 
are expressed in degrees: 


a n Cac A 
Dn = (195 E5)a+ - (=) 


Cdc = 0.1 + 1.6. TR (35) 


It will be assumed that in the most unfavourable condition 
(a = a), the value of (F/a) will be the minimum obtained in 
section 2.1. Substituting this value in expression (8) C,/a is 
obtained. With this initially assumed value of TR = 1 plus the 
C,/a value now acquired, we see that the separating angle is 
52.24° or 0.91 rad (a difference of 3.4% with respect to the 
value calculated through expression (31)). 


(34) 


Being prudent, the most unfavourable separating angle 
regarding the vessel axis is taken (0.88 rad), therefore it is 
necessary to adjust expression (35): 


C,,=0.1+1.62 TR (36) 


Because of this, we arrive at anew more conservative 
expression that defines the law for lift force per unit of flap 
rudder angle (34), which will be evaluated later. 


C,/a = 0.04+ 2.61-10+-a (37) 


Calculation of the rudder forces and profile type 


In this section, the forces per angle turn of the rudder will be 
evaluated to determine whether the minimum force demanded 
is achieved (section 2.1). For this, force coefficients will be 
analysed, of which the normal force is: 


F 


n 

Ce, /a= a 

m/ 0.5- p: V2- A, 

Taking into account the value of (F/a),,;,, calculated in 

section 2.1, with a as 0.61 rad (35°) in equation (38) (C,,/a) 

can be obtained. 

On the other hand, in fishing vessels the transverse 

components and normal force on the rudder can be connected 
in the following form [3]. 


F, =F/C (39) 


Where C is a constant determined empirically, and is 
defined as: 


(38) 


min 


C=1-0.00286u (40) 


By substituting these two expressions in equation (38), and 
bearing in mind (37), a new equation is defined, as a function 
only of turning angle a: 


Cre =" 
Cm _a 0.04 + 2.61 : 10 “a (41) 
a C 1 — 0.00286a 


In this expression, we have already remembered the selected 
rudder (with flap and taper ratio equal to 1), because expression 
(37) is used, and it was modified for this particular case. 

Although it would be convenient to evaluate expression (37) 
by comparing it with another expression of lift force coefficient 
for rudders with flaps [9]: 


C, = 2.2620 + 0.94538 = 0.93290? + 


— 0.603982 + 0.4736a8 


Where f is the angle that forms the flap with the rudder 
axis. Setting a relation between a and f that, a/p = 1 [4], for 
355 one determines that C, is 1.91 and C, is 1.95. That is to 
say, they achieve similar values; therefore expression (37) can 
be given as valid. 

So for 35° (F,/@) min = 1545.31 kN/rad (see section 2.1) 
according to equation (38). 


(C/O) min = 0.0306 
And, alternatively, with respect to equation (41) for 35s: 
C/A = 0.054 > (Cy /G) min = 0.0306 


Like this, the manoeuvrability requirements are met, in 
accordance with the normal forces and that of lift required. 

To test the rudder evolution at different angles (Fig. 8) 
expressions (8), (30) and (41) are used. The resistance coefficients 
for drag and the resulting force can be calculated by: 


(42) 


min 


POLISH MARITIME RESEARCH, No 3/2012 31 


Ca =C: cosa + C,,- sina (43) 
0.5 

HD)  # 
a a a 
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Rudder angle (degrees) 


Fig. 8. Evolution of dimensionless force coefficients 
on the rudder for different angles 


One can see in Fig. 9 that, for whatever the angle of rudder 
turn, it complies with P > 0.275. 


Rudder angle (degrees) 
Fig. 9. Number P of Norrbin for different angles 


The main characteristic that defines the type of profile is 
the relation t/c (33), as it will affect the minimum resistance 
and the stall angle. The profile characteristics are those which 
have greatest influence on cavitation. 

For the study of the flap rudder, a symmetrical NACA 0030 
profile was used (Table 2). In this profile we shall examine 
whether cavitation exists in each section, for each angle of 
action according to expressions (12) and (13). 

Firstly, the pressure coefficient (C,) has been obtained 
from expression (15), taking into account the profile surface 
perimeter (I), where pressure is expressed as: 


1= Í (1 +y’(xy)o5 dx (45) 


At the same time, the profile thickness referred to as axis 
“y” is a function of the position (x) considered on the chord 
(c). Consequently, for the NACA 0030 profile [1], y (x) is 


a percentage of c: 


tmax 
=+ , 
y= +959 
7 0.26690x%5 — 0.12600x + Jo 
—0.35160x? + 0.28430x3 — 0.10150x* 


The chord point considered for this study is the pressure 
centre for each angle of attack of the rudder. The evolution of 
C, with angle of turn is shown in Fig. 10. 

The cavitation coefficient (o) was calculated for different 
angles of attack of the rudder (a), and from different sections of 
the rudder (Table 3). The calculation of P, in each one of them 
(defined by their hg values) has been found by considering the 
draft (T) and the project trim for the vessel base (see table 1). 
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Tab. 3. Cavitation number for different transversal sections of the rudder 


Sections 


1 | The nearest to the heel pintle 165161 


Coincident with the lowest 
propeller tip 


Coincident with the 
propeller axis 


161535 


139938 


Coincident with the top 
propeller tip 


118341 


The nearest to the rudder 


107904 
root 


For every section considered and the range of rudder angles 
analysed, the risk ofexcessive cavitation is not noticeable (Fig. 
10), therefore it seems correct to assume the initially selected 
NACA0030 profile. 


3.5 


0 10 20 30 40 50 60 
a(grados) 


Fig. 10. Pressure coefficient and cavitation numbers 
by section for different angles 


Selection of rudder balancing 


The compensation (X,) is the relation between the blade 
surface to the bow axis for turning, and the total surface area. 
This parameter has a fundamental importance for reducing 
torque moments and therefore on the stock diameter. Its 
principal limitation is that the centre of pressure, considering 
that this moves with the angles of turn, always remains at the 
stern. 

The evolution of the transverse pressure centre position 
(CP,), measured as a percentage of the chord from the bow 
(see Fig. 7), for distinct attack angles has been calculated with 
expressions (36) and (41) together with [9]: 


Cmc 
CP, = | 0.25-—4 a 
fn 
Cat 
-4 = 0254 (48) 
Q 
1.11((A7+4)°°)+2 
57.3(1 +5) 
Cac Qa 
A 57.22 


Substituting the values for the base vessel (Table 1) the 
following expression is reached (with a in degrees): 


3 + 107? — 1.49 - 107ta 
0.04 + 2.61 : 1074ta 
1 — 0.00286a 


In addition, the values achieved for TR = 1 (see section 
2) should be displaced by 1% to the stern [9] (Fig. 8). The 
longitudinal centre of pressure is expressed as: 


085 o 
cr, (aa ) 


CP, = 0.25 — (49) 


(50) 


—+-CPc(m) 


0 10 20 30 40 50 60 70 
a (degrees) 


Fig. 11. Distance from the transverse pressure centre with the rudder angles 


Figure 11 shows that for angles close to zero, the distance 
from the transverse pressure centre (CP,) is 0.54 m. Therefore, 
X, = 0.54 m will be used as this guarantees positive torque 
moments. The compensation factor for the value of X, is 18%, 
and the reduction of torque respective to a compensation of 0% 
is shown in Fig. 12. 

The torque moment supported by the rudder at 35° (the 
bending moment) exceeds the minimum values demanded by 
DNV [4] for balanced rudders, as well as the recommendations 
on maximum compensable area (23%) and the maximum 
compensable chord (35%). 


4000 


3500 +— == X1=0,54 


~ X1=0 


a (degrees) 


Fig. 12. Torque moments for zero compensation and of 18% (X, = 0.54) 
for different rudder angles 


Calculation of the stock diameter 
For this calculation the following equation will be used [4]: 


1 1 
4 (MAY? Ma3 
= age fae | t (51) 
da E (=) ) (=) (mm) 


This expression requires that the torque moments (M,) are 
calculated along with the bending moments (M,), for a rudder 
angle of 35%. For these calculations expressions (50), (47), 
(38) and (44) are used. 


M, = F,((CP, < c) -X,) 
M; = F(CP/7) 


(52) 
(53) 


Finally, the diameter obtained d = 334.46 mm meets to the 
CS requirements [4]. 


Selection of rudder flap 


From the trials taken on symmetric profiles NACA for the 
rudders, with a ratio aspect, A = 2, and Re = 0.125-106[9], it 
could be concluded that the lift forces were almost double for 
rudders with a flap ratio of 0.25, and a relation of flap angle/ 
rudder angle = 2. 

A fundamental parameter, that affects the effectiveness 
of arudder flap, is the compensation factor. In section 2.5, 
a compensation of 18% was defined, and for this, according 
to experimental results mentioned, the greatest relation C,/C,, 
is given for the relation flap chord/total chord = 0.4. It is this 
value that is taken in order to proceed with the calculations. 


RESULTS 


As aconsequence of the steps indicated in the previous 
sections, an optimised rudder model, which is actively 
supported and adapted to the vessel base, has been obtained. 
This model has a NACA 0030 profile, a flap chord of 40% and 
squared edges. The rudder characteristics are summarised in 
Table 4 and its relative position is shown in Fig. 13. 


Tab. 4. Operating and geometric characteristics of the rudder 


fem] 2 [hm] ve [xt [TR] 2 | Re 


3 |2 |é | 03 [18% 1 | o [310] 


The rudder defined has been designed to comply with the 
forces per angle unit and moments, which allow the certain 
existing operational requirements of the regulations to be 
achieved, as well as those recommended by the experimental 
test results obtained (section 2). 
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Fig. 13. Dimensions and relative position 
of the rudder for tuna purse seiners (mm) 
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However, other demands are not guaranteed, such as the 


SOLAS Chapter II-1, Rule 29 extreme, or the aptitude in order 
to correct the yaw in respect to OMI Circular MSC.1053. For 
that it will be necessary to model the rudder defined here in 
a hydrodynamic design tool, which allows the model used to 
be optimised with evolutionary algorithms. 


DISCUSSION 


This paper has introduced a method of definition for the 
principal characteristics of a high-performance rudder and 
its relative position in the boat, in order to optimise the 
manoeuvrability of the fishing vessel. 

The method is initiated with the definition of the minimum 
normal, and lift forces per rudder unit angle, which must be 
exerted on the rudder to assure the “turning ability”, ‘course 
keeping’ and “yaw checking”. Integrating the results obtained 
from trials with fishing vessels and rudders have defined 
expressions that relate the rudder geometry and its relative 
position in relation to the propeller and said forces. 

Finally, the method is applied to a representative vessel with 
high manoeuvrability demands: a tuna purse seiner. The 
results obtained allow not only knowing all the parameters 
that define the rudder, but also its behaviour in various 
sections and for different angles. Also, the rudder has been 
checked for compliance to the Classification Societies’ 
requirements. 

Even though the rudder defined assures certain manoeuvrability 
conditions in the vessel, other regulative demands must be 
checked. For that, the model obtained and its respective 
evolution can be anticipated by the application of 
evolutionary algorithms which will be tested with CFDs. 
The results acquired for the base vessel show a remarkably 
high Re value reached at service speed. This, together with 
the need for manoeuvrability at large rudder angles, leads to 
a high risk of reaching early separating angles. For this, it is 
proposed that subsequent simulations of the behaviour of these 
models consider the ‘flow straightening’ effect, and that of 
cavitation along the profile surface where the laws of pressure 
show lower values than at all other respective points. 


NOMENCLATURE 
a — hysteresis loop width (m) 
A, | — rudder area (m?) 


A, — advance (m). The distance travelled by the centre of the 
vessel in the direction of the original course from the 
initiation of the turning manoeuvring until that the course 
has been modified by 90s 

B — breadth (m) 

c — average rudder chord (m) 

C — empirical coefficient which relates transversal and normal 
rudder force components 

Cb — block coefficient 

Cdc — transversal advance resistance coefficient 

Ca  — dimensionless rudder drag force coefficient 

Ca  — dimensionless rudder normal force coefficient 

Cs  — dimensionless resultant force coefficient 

Ca  — dimensionless lift force coefficient 

C,  — dimensionless pressure coefficient 

CP, — vertical distance to pressure centre (m) 

CP. — transversal distance to pressure centre (m) 

C,,./4 — torque coefficient on the first quarter of the chord 

C,  — lift coefficient for flap rudders 

d — propeller diameter (m) 

da — stock diameter (mm) 

D — turning circle diameter (m) 
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tactical diameter (m) 

course deviation (m) 

normal force to the axial plane of the rudder (N) 

lift force perpendicular to the inflow direction on the 
rudder (N) 

resulting force (N) 

material factor (see DNV rules, Pt.3 Ch.3, Sec.2) 
gravity (9.8 m/s?) 

rudder span (m) 

distance to the rudder section considered from the floating 
line (m) 

propeller advance ratio. 

coefficients of hydrodynamic absorption 

transversal radius of inertia of the mass of the vessel (in 
fishing vessels 0.24—0.26) 

inertial radius for the dragged water of the vessel (0.2 for 
the base vessel [3]) 

propeller thrust coefficient (0.18 for the base vessel with 
0.55 of blade area ratio and P/D = 1). 

length of the profile surface where the lift is acting (m) 
total length of the vessel (m) 

length between perpendiculars (m) 

displacement of the vessel (T) 

torque moment (Nm) 

root bending moment (Nm) 

propeller rate of revolutions (rps) 

Norrbin P number 

free flow pressure (Pa) 

atmospheric pressure (Pa). Pat = 101325 Pa. 

local pressure on the sections (Pa) 

local vapour pressure (Pa). P, = 1706 Pa for an average 
temperature of the sea 15sC 

Reynolds Number 

rudder section thickness (m) 

average design draft (m) 

aft draft related to the design draft (m) 

fore draft related to the design draft (m) 

taper ratio 

service speed of the vessel (m/s) 

inflow speed on the rudder (m/s) 

separation between the fore edge of the rudder and the 
propeller plane (m) 

rudder compensation (%) 

point on the chord profile 

longitudinal distance from the midship of the vessel to the 
centre of buoyancy (m) 

distance from the propeller axis to the longitudinal axis of 
the rudder (m) 

profile thickness law relative to the chord 

distance from the top tip of the propeller blade to the 
lower tip of the rudder blade (m). 


angle of attack (rad except where another unit is 
mentioned) 

angle between the axis vessel and the hydrodynamic 
inflow direction on the rudder for lift zero (rad except 
where another unit is mentioned) 

effective rudder angle (rad except where another unit is 
mentioned) 

stall angle with respect to the vessel axis 

flow straightening factor 

sea water density 1.02510 3(Kg/m3) 

kinematic viscosity of water (10° m?/s for an average 
temperature of the sea of 20sC and normal atmospheric 
conditions) 

rudder turning angle. The angle that forms the flap with 
the rudder axis. 

drift angle 

aspect ratio 

cavitation number 
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Effect of the working liquid compressibility 
on the picture of volumetric and mechanical 
losses in a high pressure displacement pump 
used in a hydrostatic drive 
Part Il 
Mechanical losses in a pump 
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ABSTRACT 


Working liquid compressibility may considerably change the values and proportions of 

coefficients of the volumetric and mechanical energy losses in the displacement pump used 

in a hydrostatic drive system. This effect can be particularly seen in the operation under 

high pressure and also in the system, where aeration of the working liquid can occur. In 

the Part II the mathematical model is presented of the torque of mechanical losses in the 

pump and its laboratory verification. Conclusions are drawn regarding the effect of working 
liquid compressibility on the mechanical and volumetric losses in the pump. 


Keywords: hydrostatic drive system; pressure displacement pump; energy losses; 
energy efficiency; energy saving; liquid compressibility; Sankey diagram 


8. MODEL OF THE MECHANICAL LOSSES 
IN THE PUMP ,, WORKING CHAMBERS - 
SHAFT” ASSEMBLY 


The pump shaft torque M, (required by the pump of its 
driving motor) must be greater than the torque M,, indicated 
in the pump working chambers because of the necessity of 
balancing also the torque M,,, of mechanical losses in the 
„working chambers - shaft” assembly. The assembly forms 
the working chambers and changes their capacity and also 
connects the working chambers with the shaft. Therefore, the 
torque M, required on the pump shaft is a sum of the torque 
M indicated in the working chambers and the torque M,,,, of 
mechanical losses in the pump ,,working chambers - shaft” 
assembly [19]: 


Mp = Mo: + Mpm (8.1) 


Torque M,,, of mechanical losses in the pump with variable 
capacity qp,, per one shaft revolution is, at the maximum value 
Of poy 1-€. poy = Gp; (with bp = qp,/qp, =1), equal to the torque 
of mechanical losses in that pump working as a pump with 
constant capacity q», per one shaft revolution. The theoretical 
and mathematical models describing the torque M,,, of 
mechanical losses in the pump with variable capacity dp, per 
one shaft revolution may be based on models of M,,,, describing 
the torque of mechanical losses in the pump with constant 
capacity qp per one shaft revolution (with b, = 1). Considering 
the models describing the torque of pump mechanical losses, 
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we assume, that, in the hydrostatic transmission, pump is driven 
with practically constant rotation speed n, and the decrease 
of shaft speed (decrease of the pump driving motor speed as 
a result of the increase of torque M, loading the motor shaft) to 
a value np < npo (Mp, — rotational speed of unloaded pump driving 
motor) is negligible from the point of view of the impact of 
speed n, on the value of torque M,,, of mechanical losses. 

Torque M,,, of mechanical losses in the pump is mainly 
an effect of friction forces between elements of the ,,working 
chambers - shaft” assembly, depending, among others, 
on the torque M, indicated in the working chambers — 
Mp; = Orgy APpi/211 = bp qr; App/2II. 

Friction forces between elements of the ,,working 
chambers - shaft” assembly are, to some extent, also an 
effect of the load on those elements of the inertia forces 
from rotational and reciprocating motion and depend 
on the pump capacity q»,, per one shaft revolution (bp 
coefficient). 

Particularly, in piston (axial or radial) pumps with 
casing (crankcase) filled with liquid, friction forces also 
occur between elements of the „working chambers - 
shaft” assembly and the liquid and depend on the liquid 
viscosity v. 

The value of torque M PmlApp;.bpv, Of mechanical losses in 
the pump ,,working chambers - shaft” assembly, loaded with 
indicated increase App, of pressure in the working chambers, 
in the pump operating at the capacity qp,, = b, qp,per one shaft 
revolution and discharging the working liquid with (constant) 


reference viscosity v,, can be described as a sum of torque 
MopmlAppiz0.bp.v, OF mechanical losses in the unloaded pump 
(torque of the losses when the indicated increase App, of 
pressure in the pump working chambers is equal to zero — 
App»; = 0) and increase AM Pelli PENs of torque of mechanical 
losses, the increase being an effect of loading the assembly 
structure elements with torque M,, indicated in the pump 
working chambers (torque M,; appearing when the indicated 
increase App; of pressure in the pump working chambers is 
greater than zero — Ap,;> 0): 


+AM (8.2) 


Merit = eal App: =0,bp,v, Pm|Appi.bp.Vn 


Torque M}; indicated in the pump working chambers is 
proportional to the increase Ap, of pressure in the chambers 
and to the active volume of the chambers created during one 
pump shaft revolution, which is equal to the theoretical capacity 
dp, (Vp) per one shaft revolution in a pump with constant 
capacity per one shaft revolution or to the geometrical capacity 
Ape = bp qr; per one shaft revolution in a pump with variable 
capacity per one shaft revolution. 

Therefore, the ,,working chambers - shaft” 
structure elements are loaded: 

— in a pump with theoretical (constant) capacity q», per one 
shaft revolution — with indicated torque: 


M is — APtAPpi Ap Pi 
211 


— ina pump with geometrical (variable) capacity qp,, per one 
shaft revolution — with indicated torque: 


ApevAPpi Z bp dp, APp; 
211 211 


which, in effect, causes a differentiated intensity of the 
increase AM of the torque of mechanical losses, 


assembly 


Mpi = 


Pm|Appi.bp.vn 
determined, with different values of coefficient b, = qp,./dp,. 
as a function of the increase Ap», of pressure in the pump 
working chambers. 


App; 
Mem bo À 
y 


Np=cte, Gpgy=0 (bp=0),qpgv (bp), apg = apr (bp =1), Vmin, I, Ymax 


In the theoretical and mathematical models describing 
the torque M BJA Dei be Na of mechanical losses a hypothesis is 
assumed, that the increase AMpalappisbp.vn of the torque of 
mechanical losses in the pump is proportional to the torque 
M»; indicated in its working chambers. 

The impact of inertia forces of the „working chambers 
- shaft” assembly elements, performing the rotational and 
reciprocating motion in the pump, on the torque M,,, of 
mechanical losses can be presented, under the assumption 
that rotational speed n, of the pump driving motor changes 
only in a small range, as a function of capacity qp,, (bp 
coefficient) per one shaft revolution of a variable capacity 
pump. Inertia forces do not depend on the value of increase 
AP»; of pressure in the working chambers, therefore their 
impact on the torque M,,, of mechanical losses in the 
pump may be included in the evaluation of the torque 
M pmlapri=0, bp,v, Of mechanical losses determined at the 
increase App; = 

The impact of the friction forces between the ,,working 
chambers - shaft” assembly elements and the liquid on 
the torque M,,, of mechanical losses in the pump can be 
presented, under the assumption that speed n, changes in 
a small range, as a relation of M,,, to the viscosity v and to 
the capacity qp,, (bp coefficient) per one shaft revolution 
(Fig. 8.1). 

It is assumed, that the impact of liquid viscosity v on 
the friction forces between the ,,working chambers - 
shaft” elements and the liquid in the piston pump casing 
(crankcase), and in effect on the torque M,,, of mechanical 
losses in the pump, can be evaluated at one level of the 
increase App; of pressure indicated in the working chambers, 
e.g. at the increase Ap,, = 0. This assumption is connected 
with a simplification assuming that there is no significant 
impact of the increase App; of pressure on the liquid viscosity 
v and with assuming in the model describing the torque Mpm of 
mechanical losses the liquid viscosity v determined in the pump 
inlet conduit (at pressure pp, equal to zero (at liquid absolute 
pressure equal to atmospheric pressure)). 


Mempo 
Vmin 
APpi= Ph 
Mpm Dp Tj 
M. á be 
1 2... 
0 Pn AP 


Fig. 8.1. Torque Mpy\y,,,, bpv 
and with variable capacity qp, = 


of mechanical losses in a pump (especially in a piston (axial or radial) pump with crankcase filled with liquid) 
bpp, per one shaft revolution, as a function of the indicated increase Ap, of pressure in the pump working chambers 


— graphical interpretation of theoretical model (9); capacity q», per one shaft revolution (coefficient bp of pump capacity): 


arge 7 0 (bp = 0), UPev (bp), dpe — 


4p, (bp = 1); liquid viscosity V ni, 


, Y, and v 


max 
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The impact of inertia forces of structure elements performing 
the rotational or reciprocating motion in the pump and also the 
impact of liquid viscosity v on the torque Mpm of mechanical 
losses in the pump is then described in the model of the torque 
M Pm|Appi=0, bp, v Of those losses in an unloaded pump (at App; = 0) 
supplied with working liquid of changing viscosity v. 

The proposed theoretical models, describing the torque 

Pm|App;=0, bp, v Of mechanical losses in an unloaded pump (at 
the indicated increase Ap,; = 0 of pressure in the working 
chambers) and at changing working liquid viscosity v, have 
the form: 

e inthe pump with theoretical (constant) capacity qp, (bp = 1) 
per one shaft revolution: 


aym 
v 
M PmjAppi=0,bp=1,v = M pr] App; =0,bp=l.v, 0 (8.3) 
n 


e inthe pump with geometrical (variable) capacity qp,,, (qpgv 
= bp qp,) per one shaft revolution: 


Pm|App;=0, bp, v = 


n 
Vi 


y yom (8.4) 
=(M pp] App;=0,bp=0, v, PAM pm|app,=0.bp. v) ci 
n 


where: 


AMpm|app;=0, bps Yn M Pmjap p¡=0, bp, vt 


=M pmjapp=0, bp=0, va (8.5) 


= (Mo mlappi=0, bp=l, vn MbmjAppi=0,bp=0, Vi bp 


Exponent a,,, in expressions (8.3) and (8.4) describes the 
impact of the ratio v/v, of working liquid v to reference 
viscosity v, = 35mms"! on the value of torque of mechanical 
losses especially in a piston displacement machine with 
liquid filling the casing (crankcase) (in the pump and in the 
hydraulic motor). 

The increase AMpni| Appi bp,v Of the torque of mechanical 
losses in the pump, due to the load of the assembly elements 
with the indicated torque M,, resulting from the indicated 
increase App; of pressure in the pump working chambers, 
is independent of the inertia forces of elements performing 
the rotational or reciprocating motion in the pump. It is also 
practically independent of the working liquid viscosity v; 
therefore, it may be determined at one viscosity value, e.g. 
at the liquid reference viscosity v, (Fig. 8.1). 

In the mathematical models describing the torque M,,, 
of mechanical losses in the pump, coefficients k; of losses 
are used relating (comparing) the components describing 
the torque M,,, of losses in theoretical models to the pump 
theoretical torque M,,. The pump theoretical torque Mp, is 
also a reference value used in the description of the torque Mp; 
indicated in the pump working chambers: 

— theoretical torque: 


_ArtPn 
Pt 


of the pump, with theoretical (constant) capacity qp, per one 
shaft revolution (b, =1), is determined with the increase 
Ap, of pressure in the pump equal to the system nominal 
pressure p, — Ap» = p,, and with the assumption that there 
are no pressure and mechanical losses in the pump, 
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— indicated torque: 


= dpc AP pi — APiPn AP pj _ AP pi 
211 


Pt 

211 py Pn 
in working chambers of the pump with theoretical (constant) 
capacity q», per one shaft revolution (b, = 1) is determined 
with the indicated increase Ap,, of pressure in the working 
chambers, 

— indicated torque: 

Miss = Poy AP pi 
i” amo 

211 


Mpi 


_ bpd pAP pi _ 
211 


APPa AP pi AP pi 
== b, —“ =M,,b 
mu > i pt Op 


in working chambers of the pump with geometrical 
(variable) capacity qp,, = bp q», per one shaft revolution is 
determined with the indicated increase Ap,, of pressure in 
the working chambers. 


The theoretical and mathematical models describe the 
torque M,,, of mechanical losses in the pump with theoretical 
(constant) capacity qp per one shaft revolution or with 
geometrical (variable) capacity qp,, = bp qp, per one shaft 
revolution: 

— pt= Aplapp;=0.pp1i=0.bp=1.v, 18 theoretical capacity per one 
shaft revolution of the pump with constant capacity per 
one revolution (bp = 1) determined at App; = 0, ppi; = 0 and 
V„ Which is equal to the working volume of the working 
chambers created during one shaft revolution, 

— drev = bp Gp, is geometrical capacity per one shaft revolution 
of the pump with variable capacity per one revolution at 
App; = 0, Ppi; = 9 and v,, which is equal to the working 
volume of the working chambers created during one shaft 
revolution. Capacity qp,, per one shaft revolution changes 
in the 0 < dp, < qr; range and coefficient bp = qp,,/qp, of the 
pump capacity changes in the 0< b, < 1 range. 


The proposed mathematical models describing the torque 
M,,, of mechanical losses in the pump, related to theoretical 
models of the torque of mechanical losses, take the form: 

e in apump with theoretical (constant) capacity qp, per one 

shaft revolution (b, = 1): 


Va 


oe App; 
Mahony Kaa Maf + ky Mp, a = 


n 


a 
= TES +k, Pe: Mp = (56) 
Va Ph 
avin A 
= KS) +ky, Poi | Hee 
Va y 211 
where: 
k, = MbmjAppi-0,bp-1 v, = Melo -0 brv (8.7) 
a My AprPn 
211 


AM 


Pm|Appi.bp=L vn _ Pm|Appi-bp=l.vn _ 


: Mp; Apr AP pi 
211 
e M PmjAppi, Dp=L VA —M pnjap p;=0,bp=1, v, _ 
dp: APpi 
211 (8.8) 
E M pal App, =p, bp=Lv, T M PmjApp;=0,bp=1, Vn _ 
4rtPn. 
211 
M 


Pm] App; =Pn bp =l, Vn M pmjAppi=0.bp=1, Vn 
Mp: 


e ina pump with geometrical (variable) capacity pa, 
(drev = bp qp,) per one shaft revolution): 


a ym 
v 
M PmjAppi.bp.v = (k41. +k4.12 bp)Mpi a + 


APpi = 


+k4.Mp;,bp 
n (8.9) 


Avm 
v Apri 
=[(Ky ¡1 tky 12 bof | +k,» bp 7 Mp = 


n n 


where: 


aym 
v Appi |qptPn 
= (Kg 11 +K412Dp) | +Kg42 bp — |= 
n n 211 
k E Mpmlappi-0,bp=0, Vn Mo mnlappi-0, bp—0,Vp (8.10) 
Aq i ~r O ` 
Mp arPr 
211 
k _ M pal App; =0,bp=L, Vn Milo, bp=0.Vp 
412 7 M = 
j (8.11) 
MBmjAppi=0.bp=1, v, —Mpmjapp;=0, bp=0,v, 
APP in 
211 
— AMomn|Appi.bp. vn 2 AMomlappi.bp.Vn _ 
A TO Aa 
Mp, bp dpAPpi 
21 
_ AM PmjAppibp=1, v, _ 
ApAPpi 
211 (8.12) 
Mo mlape:=Dasbp=l.Vy = Monlspp.=0,bp=l,¥, _ 
QptPn 
211 


= M bmjAppi-pn.bp-1, vn = Mpmjappi=0,bp-1, v, 


My 


Commentary: 

— The sum (k,,,+k,,,) of coefficients used in mathematical 
model (8.9) describing the torque M,,, of mechanical losses 
in the pump with geometrical (variable) capacity Qp, (Gpey 
= bp q») per one shaft revolution is equal to coefficient k, , 
used in the mathematical model (8.6) describing the torque 
Mpm of mechanical losses in that pump working as a pump 
with theoretical (constant) capacity per one shaft revolution: 
Kara + ky = ky. 

— Coefficient k,,used in mathematical model (8.9) describing 
the torque M,,, of mechanical losses in the pump with 
geometrical (variable) capacity qp,, (Upa = bp qp) per 
one shaft revolution is equal to coefficient k,, used in 
the mathematical model (8.6) describing the torque Mpm 
of mechanical losses in that pump working as a pump 
with theoretical (constant) capacity qp, per one shaft 
revolution. 


9. RESULTS OF INVESTIGATION OF THE 
TORQUE M,,, OF MECHANICAL LOSSES 
IN THE MEDIUM PRESSURE AXIAL 
PISTON VARIABLE CAPACITY PT0Z2-25 
TYPE PUMP USED IN A HYDROSTATIC 
TRANSMISSION SYSTEM 


Michał Czyński [21] performed the investigation of the 
PT0Z2-25 medium pressure pump as laboratory verification of 
the mathematical model of the hydrostatic transmission energy 
efficiency. The pump was loaded by a hydraulic motor and not 
by an overflow valve (as it was in the case of the investigation 
of the high pressure axial piston HYDROMATIK A7V.58. 
DR.1.R.P.F.00 type pump). The threat of the hydraulic oil 
aeration was smaller. The system operated under the nominal 
pressure p, = 16MPa, half of the value of pressure used in the 
investigation of the A7V.58.DR.1.R.P.F.00 pump. 

The investigation lasted relatively short time, the hydraulic 
oil tank was large in relation to the pump power. Therefore, 
the effect of oil aeration and of its compressibility was 
minimized. 

The torque M,,, of mechanical losses in the „working 
chambers - shaft” assembly of the pump is determined from 
the measurement of the shaft torque M, and the torque Mp; 
indicated in the working chambers. Direct measurement of 
the torque M,, is not possible. Therefore, it is necessary to 
determine the torque Mpm of mechanical losses by an indirect 
method from the formula: 

pey APp; 
Mpm =Mp —Mp; =Mp s= = 
211 
bp qp: AP pi 
211 

Fig. 9.1 presents results of the investigation of the torque 
M, on the pump shaft as a function of the indicated increase 
App; Of pressure in the working chambers for different values 
of the coefficient b, of the pump capacity qp,, = bp qp; per one 
pump revolution. The investigation results show a very high 
(close to 1) value of the determination coefficient R2. 

Fig. 9.2 and 9.3 present the calculated values of the torque 
M,, indicated in the pump working chambers and the torque 
Mpm Of mechanical losses in the pump „working chambers 
- shaft” assembly. 

Fig. 9.4 presents the torque M pmjAppi=0,bp.v, Of mechanical 
losses in the unloaded pump as a function of the pump capacity 
coefficient bp. 


(9.1) 
=M, - 
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Fig. 9.1. Torque M, on the PTOZ2-25 pump shaft as a function of the Fig. 9.3. Torque M,,, of mechanical losses in the PTOZ2-25 pump ,, working 
indicated increase App; of pressure in the working chambers, for different chambers - shaft” assembly as a function of the indicated increase App, of 


values of the pump capacity coefficient bp at v, working liquid viscosity [ 21] pressure in the working chambers, for different values of the pump capacity 
coefficient bp at v, working liquid viscosity [21] 


bp=1.000; R?=0,999 


48 bp=0.945; R?=0.999 | 
—————— b,=0.897; R?=0.999 

44 "== b,=0.849; R?=0.999 
——————- b,=0.795; R?=0.999 

40 —- ———— bp=0.737; R?=0,999 


bp=0.690; R*=0,999 
b,=0.593; R*=0.999 
b,=0.500; R?=0.999 
b,=0.397; R?=0.999 
b,=0.299; R?=0.999 


torque M,, indicated in the working chambers [Nm] 
8 


torque Momsp,,,-0 Of mechanical losses in the unload pump [Nm] 
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Fig. 9.2. Torque M,, indicated in the PTOZ2-25 pump working chambers Fig. 9.4. Torque Mpjap,=0 yv, of mechanical losses in the unloaded PTOZ2- 
as a function of the indicated increase App, of pressure in the working 25 pump „working chambers - shaft” assembly as a function of the pump 
chambers, for different values of the pump capacity coefficient bp at v, capacity coefficient b, at v, working liquid viscosity [21] 


working liquid viscosity [21] 
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Fig. 9.5 presents the values of increase AM pmjappi,bp.v, OF 
torque of mechanical losses calculated from the investigation 
results. The increase AM pmlappi.bp.¥n is proportional to the 
increase App, of pressure in the working chambers and to 
the instantaneous value qp,, = bp qp, of the pump geometrical 
working volume. 


2.2 —————— b,=1.000; R*=0.983 — 
b,=0.945; R?=0.984 
20 b,=0.849; R?=0.977 | | | ? 
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1.6 — —————— b,=0.299; R?=0.942 + TIAS 


increase AM,,, of torque of mechanical losses in the pump [Nm] 
> 


0 2 4 6 8 10 12 14 16 
indicated increase Ap,, of pressure in the working chambers [MPa] 


Fig. 9.5. Increase AMpjap,,, r.v, Of the torque of mechanical losses in the 
PTOZ2-25 pump ,, working chambers - shaft” assembly as a function of the 
indicated increase App; of pressure in the working chambers, for different 
values of the pump capacity coefficient bp, at v, working liquid viscosity [21] 
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Fig. 10.2. Picture of the torque Mp, 


Results of the investigation of the torque Mp, of mechanical 
losses in the PTOZ2 - 25 pump ,,working chambers - shaft” 
assembly confirm the mathematical models ((8.9), (8.12)) 
of mechanical losses in the displacement pump at reference 
viscosity v, = 35mm?s"!. 

The above presented results can be achieved with 
the assumption of little hydraulic oil compressibility, i.e. 
compressibility with insignificant effect on the change 
of pump theoretical working volume q», (Vp) of the 
working chambers (or the geometrical working volume 
drev = bp q») under the influence of the change of indicated 
increase App; of pressure in the working chambers. 

The picture of the increase AMpni| Appisbp.v, Of the torque 
of mechanical losses in the pump ,,working chambers - 
shaft” assembly, complying with the mathematical model 
(equations (8.9) and (8.12), Fig. 9.5), can be a tool for 
evaluation of the coefficient k,,,, of compressibility of the 
working liquid used in the hydrostatic transmission system. 
The picture of AM Pm|App;,bp.v, different from that model is then 
a result of the liquid compressibility, which reduces the values 
pil appi =P, (OF dPevlap 9 bp Gp: Appi =P) under the influence of 
the increase App; and, in consequence, reduces the value of 
indicated torque Mp; in the working chambers and torque M, 
on the pump shaft. 

In consequence, the value of the increase AM pap; bp y, OF 
the torque of mechanical losses, calculated from expression 
(9.1), decreases as a difference between torque M, measured on 
the pump shaft and torque M,, calculated with the assumption 
of non-decreased value q»; (Vp) (Or pg, = bp Gp). Also decreases 


the value of coefficient k, , of the increase AM lb 


10. RESULTS OF INVESTIGATION OF THE 
TORQUE M,,, OF MECHANICAL LOSSES 
IN THE HIGH PRESSURE PISTON PUMP 

HYDROMATIK A7V.DR.1.R.P.F.00 TYPE 
LOADED BY AN OVERFLOW VALVE 


Jan Koralewski [22], apart from the verification investigation 
of the model of volumetric losses of the axial piston variable 
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of pressure in the working chambers, at the constant value v/v, = 1 of the ratio of working liquid viscosity v to the reference viscosity v,, at the pump capacity 
coefficient b, and at different assumed values of working liquid compressibility k| 324pa: 0.030; 0.035 (pump of the HYDROMATIK A7V.DR.1.R.P.F.00 type) [22] 
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Fig. 10.1. Picture of the torque M,,, of the pump mechanical losses in the pump ,, working chambers - shaft” assembly as a function of the indicated increase 
App; of pressure in the working chambers, at the constant value v /v, = 1 of the ratio of working liquid viscosity v to the reference viscosity v,, at different 
values of the pump capacity coefficient b, and at different assumed values of working liquid compressibility koi ¿»wpa 0.000; 0.010; 0.020; 0.025 (pump of the 
HYDROMATIK A7V.DR.I.R.PE00 type) [22] 


displacement pump of bent axis design (HYDROMATIK A7V. 
DR.1.R.P.F.00) (presented in chapter 7), performed, on the 
test stand presented in Fig. 7.2, the investigation verifying the 
models (8.9), (8.10), (8.11) and (8.12) of torque of mechanical 
losses in that pump. 

Investigations of volumetric and mechanical losses were 
carried out simultaneously. The pump was loaded with an 
overflow valve. 

As the investigation of torque Mom Appzo Of mechanical 
losses in the „working chambers - shaft” assembly of the 
unloaded pump (at App; = 0) confirmed fully the (8.9), 
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(8.10), and (8.11) models, the investigation of the increase 

PmlApp;.bp.v, Connected with the indicated increase App; of 
pressure in the pump working chambers did not confirm the 
model (8.12). 

Picture of the increase AM pinlAppi.bp.vn of the torque of 
mechanical losses determined by indirect method from 
expression (9.1) at the assumption, that the coefficient of 
working liquid compressibility is equal to zero (k,.\3:p. = 
= 0.000), presented in Fig. 10.1, differed significantly from 
the expected picture from models (8.9) and (8.12) and from 
Fig. 9.5. 


Calculations performed with the assumed coefficient of 
working liquid compressibility k,,| ,4p, > 0 showed the picture 
(Fig. 10.1, Fig. 10.2) of the increase AM closer to 
models (8.9) and (8.12) and to Fig. 9.5. 

Picture most similar to those of models (8.9) and (8.12) is 
diagram AM PmlApp;.bp.v, Obtained with the assumed coefficient of 
working liquid compressibility k,, | 3p, = 0.030 (Fig. 10.2). 

With the assumption of the coefficient k, |32mpa = 0.030, the 
obtained increase of the coefficient k,, of AM pmjapni bpv, Of the 
torque of mechanical losses was from the value k,, = 0.014 
to the value of the order of k,, = 0.045. 

At the same time, taking into account the compressibility 
of working liquid defined by the coefficient k,.| 3p, = 0.030, 
the decrease is obtained of the coefficient k, of volumetric 
losses in the pump (expression (7)) from the value k, = 0.065 
to the value of the order of k, = 0.035. 


CONCLUSIONS 


Pm|Appi.bp.vn 


1. Energy investigation of the high pressure displacement 
variable capacity pump operating in an aerated hydrostatic 
system allows to determine approximate value of the 
coefficient kilp, of working liquid compressibility by 
comparing the picture of the increase AM Pmlapp;,bp.v, Of the 
torque of mechanical losses in the ,,working chambers - 
shaft” assembly with the model of increase of that torque. 

2. The approximate determination of the value of k,; 
coefficient of working liquid compressibility allows i 
modify the coefficient k4, of the increase AM pnjappi.bp. yy OF 
the torque of mechanical losses in the pump „working 
chambers - shaft” assembly (towards its increase) and the 
value of coefficient k, of volumetric losses in the pump 
working chambers (towards its decrease). 

3. Therefore, it is possible to determine approximately the 
effect of the working liquid compressibility on the picture 
of values and proportion of the mechanical losses and 
volumetric losses in the pump. 

4. Working liquid compressibility reduces the pump theoretical 
working volume qp,| Appi“ Bn of the working chambers (or the 
geometrical working volume rev | App; =P, = bp qp,| bp; -»,) 
under the influence of the change of indicated increase App; 
of pressure in the working chambers. 

5. The simplifying assumption is adopted, that the sum of the 
values of working liquid compressibility coefficient Kelp, 
and coefficient k, of volumetric losses in the pump equals 
to the value of gocfivient k, when it includes also the effect 
of working liquid compressibility: 

6. Inthe mathematical models of pump energy efficiency and 
hydrostatic drive system efficiency and also of the system 
operation field (0 < OM < Duma 0 < Mm < Mmmax) either 
the coefficient k, (which includes also the effect of working 
liquid compressibility) or the above defined sum (k,|, + 

+ k,) must be used. ” 
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Computational model for simulation 
of lifeboat motions during its launching 
from ship in rough seas 


Pawet Dymarski, Ph.D. 
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Gdańsk University of Technology 


ABSTRACT 


This paper presents a computational model which describes motion of an object of six degrees of freedom 

(DoF), intended for simulation of spatial motion of one- or two- rope-sling lifeboat or rescue boat during 

its launching from ship in rough sea. This is a complex model which accounts for sea conditions as well 

as elasticity and damping properties of davit’s elements and mechanisms, rope and boat hull. Also, are 

presented results of example calculations for an assumed set of technical parameters of davit and boat as 
well as sea conditions. 


Key words: ship in rough seas; life saving appliances (LSA); lifeboat launching; 
lifeboat motions; lifeboat launching simulation 


INTRODUCTION 


Human life protection and rescue is one of the most 
important problems which must be faced by designers of objects 
doing duty for man, especially sea-going passenger ships. 

Life saving appliances (LSA) systems which make it possible 
to evacuate persons especially from large passenger ships, have 
evaluated beginning from the simplest open lifeboats and rafts 
thrown off into water and ending to unsinkable sheltered boats 
of high strength and fire resistance, lowered with the use of 
more and more advanced side davits. 

In Fig. 1 is presented a photo of an example contemporary 
davit together with suspended boat during launching it into 
water. 

Operation of launching boats together with embarked 
persons into water is the most hazardous phase of process 
of rescuing persons from endangered ship in rough seas. 
Boat, during its launching from a significant height in the 
neighbourhood of rolling ship’s side, often impacts against 
the side. In the case of large passenger ships the lowering 
height is significant and boat displaces close to ship side, 
that often leads to its bumping against the side. The bumps 
generate the relatively largest overloading which results from 
a change of motion of the boat and affects persons inside the 
boat. The overloads may lead to failures of boat and person 
injures and even casualties. For this reason many recognized 
research centres as well as leading producers of ship life 
saving and rescue appliances are searching for more and more 
exact methods for calculation of motion parameters of boat 
lowered from ship to water in rough seas conditions. When 
a reliable computational software is at disposal it is possible to 
perform investigations of LSA systems intended for rough seas 
conditions. This allows to improve the existing design solutions 
and test new ones of the devices in question. 


l 


Fig. 1. View of the lifeboat and the luffing-jib davit during launching test 


The problem was included into the scope of the research 
program realized, in the period of 2004-2009, within the frame 
of the European project SAFECRAFTS (Safe Abandoning of 
Ships - improvement of current Life Saving Appliances Systems) 
in which also the Faculty of Ocean Engineering and Ship 
Technology, Gdansk University of Technology, took part [2]. 
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Fig. 2. Coordinate systems: x, y, z — absolute coordinate system, Xy yy Z — local coordinate system of the boat, 
Xp Ye Zg — local coordinate system of the ship 


This paper presents results of one of the project's tasks 
realized by these authors, which concerns elaboration of 
a complex mathematical model and computational program, 
called RESBO, intended for simulation of lifeboat motion 
during its launching from ship in rough seas. 


GENERAL DESCRIPTION 
OF RESBO PROGRAM 


RESBO program serves to simulation of 3D motion of one- 
or two — suspension - point object during its launching from 
the ship in waves and its recovering as well. The program is 
prepared for simulation of the 1 - point lifeboat, 1- point fast 
rescue boat, 1- point davit for launching life raft and 2 - point 
lifeboat. 

The RESBO was written in C++ language and compiled 
under Linux system. 

Algorithm of the program takes into account: 

e movement of the ship in waves 

e operation of davit mechanisms 

e elasticity and damping properties of davit parts 

e elasticity and damping properties of object hull 

e influence of wind 

e lifeboat impact against ship’s side 

e lifeboat impact against water during launching 

e movement of the lifeboats in waves during a short period 
after releasing the hooks. 


THEORETICAL MODEL 
Assumptions 


The object (e.g. lifeboat) is modelled as a rigid body. 
However special flexible elements called “fenders” are attached 
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to one side of the lifeboat hull. The object performs motions of 

6 DoF. Elements and factors restricting (extorting) motion of 

the object are formulated as external forces dependent on_ 

e relative position between the object and a restricting 
(extorting) element 

e relative velocity 

e time 


Influence of the elasticity and damping properties of the 
object structure on its motion is formulated as a separate 
segment of the extorting forces function. 

The elements and factors restricting (extorting) motion of 
the object are the following: 

* suspension system and mechanisms: ropes, davit arm, 
absorber, brake, with regard to elasticity and damping 
properties of the hull construction in hook suspension 
point 

e  ship’s side which restricts area range of the object’s motions 
(possibility of the impact against ship’s side) 

e wind — aerodynamic reaction 

e force of gravity 

e water — hydrostatic and hydrodynamic reaction, including 
slamming 

e functioning the hydrostatic locking device 


The launching system includes: brake and the elements 
connected in series: davit arm, absorber, rope and catch. The 
elements transfer the same load, hence deformation of the 
system is a sum of their deformations. 


The object’s motion equations 


The dynamic equations of free spatial motion of the solid 
are as follows [1]: 


k 

he Px, Z Ce, Lao Py, Pa, ~ y Mos 
t=1 
k 

ly, Py, ~ Vay — Ixy) Pz, Po, E Y Moi 
i=1 


k 
la Pa, 7 Us, Ayo Pe, Po, = > Ma 
t=1 


Xo Vaz — coordinates of the solid’s mass centre 

Xo» Vor Zo — principal central axes of inertia (moving 
system) 

x,» Py,» Pa — Momentary angular velocities of the solid 
rotating around principal central axes of 
inertia 

JxJy Ja ~ principal central moments of inertia of the 
boat 

Fo Fy Fæ — vectors of „i” force 

M, M, , M,,— vectors of moment of ,,1” force relative to 


central axes of inertia 


Forces 


The forces which influence the object F, depend on: 

e the position of the gravity centre of the object, x,, and the 
vector of the angular position, @, (position of the object’s 
local coordinate system determined by the principal central 
axes of inertia: X, Yo Zo) . 

e the vector of linear velocity, x,, and angular velocity, € 

e the wind velocity and direction, Vy ing 

e timet 


F; = Fio že Ps Ps Vuina» t) 
Detailed description of the forces function F; follows from 
reaction of the elements on the object. 


Characteristics of the system’s elements 
Davit 


The davit is modelled as a linear springy element (without 
restriction of linear range). The stiffness of the davit, EI, is an 
input parameter of the algorithm. 


Absorber 


The absorber is modelled as a linear springy element with 
restrictions for the deformation: Al, and Al, As assumed, the 
absorber starts working when the force in the rope, Fr, reaches 
its nominal value Frn, and stops when it obtains the value of 
1,5 Frn. The absorber has immanent friction — 1f velocity of 
the deformation is not zero, then the friction force Ta which 
inhibits progress of the process, appears. Module of the force 
has constant value, independent of the force transferred by the 
absorber. The absorber is defined by the following parameters: 
Alo, Fr, Ta, vo. 

The parameter v, - velocity of friction waning (below this 
value of deformation velocity, friction linearly decreases to 
zero) is added only to obtain numerical solution. It does not 
follow from physical model of the system. 


AL 


ALa 


Frn 1.5 x Fm Fr 


Fig. 3. Assumed static characteristic of the absorber 


Steel rope 


The rope is modelled as a linear springy element with the 
minimum restriction for the transferred force: Q, = 0. It means 
that the force in the rope cannot be negative, but process of 
the rope length reduction can be still in progress. The rope is 
defined by the following parameters: its length 1, cross-section 
area A, Young modulus E. 


“Catch” of the rope to the object 


This is an imaginary element added for taking into account 
elasticity and damping properties of the object’s hull structure. 
When we simulate a system without absorber, function of 
the catch can be especially noticeable. The catch is modelled 
as a linear springy element. Its damping is modelled as an 
immanent friction, linearly dependent of transferred force. 

The catch characteristic is defined by the following 
parameters: its stiffness I, friction coefficient f, velocity of 
friction waning v, 


Fender 


The fender is modelled as a linear springy element with 
a minimum restriction for the transferred force. It cannot work 
on stretching. It starts working in the moment of contact with 
ship’s side. Its damping properties are: immanent friction which 
is in opposite direction to the perpendicular vector of the fender 
velocity relative to ship side surface. Damping force depends 
linearly on fender pressure force. Tangent friction force of the 
fender on the ship side surface is modelled analogously. 

The fender characteristic is defined by the following 
parameters: its stiffness I, perpendicular friction coefficient f,, 
tangent friction coefficient f, velocity of friction waning vo. 


Centrifugal brake 

The centrifugal brake controls unreeling velocity of the 
ropes during launching the boat. The velocity is linearly 
dependent of force transferred by the rope. The brake’s 
characteristic is presented in Fig. 4. 

Hydrostatic locking device 

The hydrostatic locking device is used for releasing the hook 

and disconnecting lifeboat’s steel ropes during launching. The 


device is defined by its position level. 


Aerodynamic reaction 


Motion of the object in the air causes aerodynamic reaction. 
This force can be described as follows: 
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Vr 


Vrn 
Vro 


Fro Ern F 


Fig. 4. Assumed characteristic of the centrifugal brake, where: 

Vro — the minimum value of the velocity Vr of the ropes which triggers the 
centrifugal brake and begins to control of unreeling velocity of the ropes 
during launching the boat, Fro, Frn — the minimum and nominal values 

of the load in the ropes during controlled launching the boat 


olV, IVS 

Ris c, Pavas alva 

where 

CA- air resistance coefficient 

p — air density 

VA-  object's velocity relative to the air, including wind 


windage area 
Hydrostatic reaction 


The hydrostatic reaction is connected with hydrostatic 
pressure exerted on lifeboat hull’s wetted area. The reaction 
which acts on the area element dS in the direction determined 
by the unit vector n, is calculated as follows: 


dF, =— py g h n dS 


where: 

P. — water density 

g — gravitational acceleration 

h — immersion depth of the area element dS 
n — unit vector normal to the area element dS 


The total hydrostatic reaction is calculated as follows: 


== [Pagh nds 
S 


Hydrodynamic reaction 


Lifeboat movement in the water induces hydrodynamic 
reaction which is roughly equal to sum of frictional resistance 
and pressure resistance. The force vector of the frictional 
resistance is calculated as follows: 


Re = Cp” Py fiolo edt dS 
S 


where: 

C, — frictional resistance coefficient 

v, — tangential component of the velocity vector of the area 
element dS, in relation to the water 


t — unit vector tangent to the area element dS, defined as 
follows: 
vy —(v, n)n 
D == RÁ 
t 

lv w —@y mal 
where: 
n — unit vector normal to the area element dS. 
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The force vector of the pressure resistance is calculated 
as follows: 


R, = C,” Pw fo, - ni(v, - n)n dS 
S 


where: 
C- 


a pressure resistance coefficient. 


Free surface reaction (slamming) 


The pressure caused by hull impact against free water 
surface is calculated by using the von Karman formula: 


_ Pudo  Tetg(a) 
Pslam = > i ae 
—2W 
where 
v, — first impact velocity [m/s] 
a — angle of inclination of the lifeboat bottom surfaces 


[deg] 
Py — water density [kg/m3] 
g — gravitational acceleration [m/s?] 
x — breadth of the wetted part of the bottom [m] 
W — weight of the body per unit length [N/m] 


The total reaction of slamming is calculated by using the 
following integral: 


Raa "JP n dS 


where: 
n — unit vector normal to the area element dS 


NUMERICAL MODEL 


The applied calculation algorithm of the object motion is 
based on the Euler explicit method. Its essence can be described 
as follows: 

0. START. Initiation of the constant values. Reading input 
data. 

1. Number of updating: n = n+1 

. Calculation of current position of ship. 

3. Calculation of forces in the elements according to positions 
of ship and object. 

4. Summing the forces and reactions which influence the 
object: 

Gravity force 

Aerodynamic reaction 

Force in the rope 

Force in the fender 

Hydrostatic reaction 

Hydrodynamic reaction 

Force of the free water surface impact (slamming) 

5. Calculation of anew position and velocity of the object, 

resulting from forces working in the time interval At. 

Calculation of a new length of the rope 

Cage the time: t= t,+ At. 

Ift then go to STOP; if not, go to 1. 


RESBOVI PROGRAM 


mono sp 


ga 


GOSS 


n+ i~ trax, 


The RESBOVI program is the post-processor for the RESBO 
program. It serves to visualize (animate) the calculation results 
of the lifeboat (or rescue boat) motion during its launching (and 
recovering) from the ship in waves. The RESBOVI was written 
in C++ language and compiled under Linux system. 


Fig. 5. Lifeboat during launching. Freeze frames (dumps) taken from the RESBOVI animation 


In Fig. 5 are given 4 successive screen dumps of animation 
which presents calculation results based on the following 
example test data:_ 

+ Heave and roll period of ship: T = 8 s 

+ Heave amplitude: Z, = 1.0m 

¢ Drift amplitude: Y, = 0.8 m 

+ Roll amplitude: d,=1.0m 

e Distance between ship’s side and its symmetry plane (PS): 

B/2=16.1m 
e Distance between boat’s centre of gravity(CG) and ship’s 

PS: bl = 20.1 m 
e Distance between ropes suspending the boat: 1 = 8 m 
e Initial position of the suspension point: h = 20 m 
e Mass of the boat: m = 15550 kg 
e Boats moments of inertia: Jx, Jy, Jz = 36792, 80102, 

80102 kg, respectively 
e Steel rope diameter: d = 16 mm 
e Nominal lay-out velocity of rope: V,,, = 0.8 m/s 


CALCULATION CASES 


The calculation cases concern the two- suspension- point 
lifeboat in 3D movement during launching (and recovering) 
from the ship in waves. Two cases are selected to compare 
their calculation results. 

The first case applies to the system with elasticity and 
damping properties of the davit arm, rope, centrifugal brake 
and boat’s hull. For comparison purposes, a small absorber is 
added to the system which is used for fast rescue boat davit. 
The second case, without any absorber, applies to the lifeboat 
davit system. Below in the sections 6.1 and 6.2 in the successive 
figures are presented calculation results for both the considered 
systems in the form of runs of the main parameters which 
describe motion of the life boat with embarked persons during 
its launching from ship to water in rough seas. The successive 
figures present the following: 
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- Fig. 6 and 11: Displacement (linear motion) of the boat's 
centre of gravity against the Earth (in the coordinate system 
x,y,z — acc. Fig. 2). 

Gx — displacements along x- axis — horizontal parallel 
to the ship’s plane of symmetry (PS); 

Gy — displacements along y- axis — horizontal 
perpendicular to the ship’s plane of symmetry (PS); 
G, — displacements along z- axis — vertical. 

- Fig. 7 and 12: Run of change in position of the local boat- 
fixed coordinate system in the global system: 

Exo,, Exo,, Exo, — coordinates of the unit vector Exo 
which determines direction and sense of Xo — axis of 
the local boat-fixed coordinate system; 
Eyo,, Eyo,, Eyo, — coordinates of the unit vector Eyo 
which determines direction and sense of Yo — axis of 
the local boat-fixed coordinate system; 
Ezo,, Ezo,, Ezo, — coordinates of the unit vector Ezo 
which determines direction and sense of Zo — axis of 
the local boat-fixed coordinate system). 

- Fig. 8 and 13: Displacement (linear motion) of the centre 
of gravity of boat and its oscillations (angular motions) 
against ship’s side (the coordinate system Xy, Y y, Zz — acc. 
Fig. 2). 

Gy, — displacements along X, - axis —parallel to the 
ship’s plane of symmetry (PS); 

Gyg — displacements along Y, - axis —parallel to the 
ship’s plane of symmetry (PS); 

G,, — displacements along Z, - axis —perpendicular 
to the ship’s plane of symmetry (PS); 

Roll, pitch, yaw — angular displacements (motions) 
of the boat; 

- Fig. 9 and 14: Linear and angular accelerations acting onto 
the boat. 

a, a, a; components of linear acceleration vector; 
£, €, €, — components of angular acceleration around 
the axis X,Y, Z, respectively; 

- Fig. 10 and 15: Linear and angular velocities of the boat 
against the Earth. 

Vo Vp V,— components of linear velocity vector; 
0, 0, 0,— components of angular velocity vector. 


System with absorber and elasticity 
and damping properties 
Diagrams 


Absolute motions of the boat 
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Fig. 6. Linear and angular motions of the lifeboat. 
Absolute coordinate system. Diagram of the mass centre motions 
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Absolute motions of the boat 
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Fig. 7. Angular motions of the lifeboat. Absolute coordinate system. 
Unit vectors of the boats local coordinate system 


Relative motions of the boat 
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Fig. 8. Linear and angular motions of the lifeboat. Ship side 5 coordinate 
system. Diagram of the relative motions of lifeboats mass centre and the 
lifeboat rotations: roll, pitch, yaw 


Linear and angular accelerations 
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Fig. 9. Linear and angular accelerations of the lifeboat 


Linear and angular velocities 
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Fig. 10. Linear and angular velocities of the lifeboat 


System without absorber 
Diagrams 


Absolute motions of the boat 
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Fig. 11. Linear and angular motions of the lifeboat. Absolute coordinate 
system. Diagram of the mass centre motions 
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Fig. 12. Angular motions of the lifeboat. Absolute coordinate system. Unit 
vectors of the boats local coordinate system 


Relative motions of the boat 
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Fig. 13. Linear and angular motions of the lifeboat. Ship side s coordinate 


system. Diagram of the relative motions of lifeboat 5 mass centre and the 
lifeboat rotations: roll, pitch, yaw 
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Fig. 14. Linear and angular accelerations of the lifeboat 


Linear and angular velocities 
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Fig. 15. Linear and angular velocities of the lifeboat 
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DESCRIPTION OF RESULTS AND 
CONCLUSIONS 


The graphs in Fig. 6 and 11 show position of lifeboat’s mass 
centre in the absolute coordinate system (x, y, z). The graphs 
in Fig. 8 and 13 show relative, linear and angular motions 
of the lifeboat in ship side’s co-ordinate system. The graph 
G z B shows sway. During launching (lengthening the 
rope) we can notice increasing the period and amplitude 
of the motion. At about 18" second of launching process 
we observe impact of the lifeboat against ship side. 

The graphs in Fig. 7 and 10 show positions of unit vectors 
of the local lifeboat coordinate system in the absolute 
coordinate system. 

The graphs in Fig. 9 and 14 show linear and angular 
accelerations of the lifeboat. The blue ones (a_z) show 
vertical accelerations. The green ones show horizontal 
accelerations (a_y). At about 4 second we can notice start 
of lifeboat lowering, and at about 18" second - impact of 
the lifeboat against ship side. 

The graphs in Fig. 10 and 15 show linear and angular 
velocities of the lifeboat. 

The presented calculation results concern a typical 
launching operation of lifeboat for two cases which differ 
to each other with rope system only. The first case concerns 
the system equipped with an absorber and the other — that 
without absorber. In view of that the absorber works 
only under load in the rope greater than nominal one the 
courses of the functions are reliable and the differences 
between two presented cases are inconsiderable. In the 
graphs in Fig. 9 and 14 the main expected difference in 
lifeboat vertical accelerations can be observed. However 
it should be stressed that in the presented model elasticity 
of davit’s structural elements was taken into account but 
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their damping properties were neglected as to estimate them 
was difficult. If the properties are taken into account, then 
values of calculated accelerations will be lower, hence the 
difference between the acceleration courses for the two 
calculated cases will be smaller. 

The program in question is prepared for further development. 
It is possible to add elements for simulating more complex 
absorbers (e.g. dampers of longitudinal and transversal 
motions of the boat or drive and control systems for 
compensators). It could be useful in elaborating new 
systems. 
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Genetic algorithm for solving celestial 
navigation fix problems 


Ming-Cheng Tsou, Ph.D. 
National Kaohsiung Marine University, Taiwan 


ABSTRACT 


As we enter the 21% century and advance further into the information age, traditional methods for computing 
a celestial navigation fix can no longer meet the requirements of modern vessels in terms of calculation 
speed and precision. Study of precise, rapid, and convenient celestial navigation computational methods 
and the application of information technology to modern celestial navigation is especially meaningful, 
considering the current push for e-Navigation. In this work, we employ a genetic algorithm, from the field 
of artificial intelligence, due to its superior search ability that mimics the natural process of biological 
evolution. Unique encodings and genetic operators designed in this study, in combination with the fix 
principle of celestial circles of equal altitude in celestial navigation, allow the rapid and direct attainment 
of accurate optimum vessel position. Test results indicate that this method has more flexibility, and avoids 
tedious and complicated computation and graphical procedures. 


Key words: genetic algorithm; celestial navigation; intercept method; celestial fix 


INTRODUCTION 


Celestial navigation is one of the important positioning and 
navigation methods for ocean voyages due to its low cost, high 
reliability, strong independence, and inability to be detected. 
It is also one of the traditional navigational skills required of 
navigators by the International Maritime Organization (IMO). 
However, since 1994, the revolutionary Global Positioning 
System (GPS) has become widely available to the general 
public. GPS meets almost all of the requirements of seafarers, 
such as high precision, all weather functionality, automatic 
operation, and global coverage. This has resulted in the 
successive abandonment of other radio navigation systems 
(e.g. OMEGA, NNSS, DECCA, and LORAN A). The celestial 
navigation system, with its longstanding history, also faces 
a severe situation. However, while we enjoy the convenience 
that GPS has brought us, we are also exposed to the risks 
arising from our over-reliance on it. Navigation safety is 
endangered when satellite signals are cut off or interfered with, 
or when the satellite receiver system is broken; both situations 
render positioning of the vessel impossible. Even today, when 
navigation is dominated by GPS, a celestial fix still serves as 
an important backup measure. Nonetheless, as we enter the 
21st century and advance into the information age, traditional 
methods for computing a celestial navigation fix can no longer 
meet the requirements of modern vessels in terms of calculation 
speed and precision. 

At the Manila Conference hosted by the IMO in 2010, a new 
amendment regarding the Seafarers’ Training, Certification, 
and Watchkeeping (STCW) Code was issued. The amendment 


describes e-Navigation as the focus of the development of 
navigation technology, and the application of information 
technology as akey direction of development. In addition, 
while celestial navigation received sustained emphasis in the 
Manila amendment, there is an explicit encouragement for the 
usage of electronic nautical almanacs and celestial navigation 
calculation software in Section B-II/1 of Chapter 2, “Guidance 
regarding the master and the deck department”, which concerns 
training and competency requirements. This can be seen as an 
effort to lead the development of traditional celestial navigation 
gradually toward e-Navigation. 

For these reasons, the study of high precision, rapid, and 
convenient celestial navigation computational methods and 
the application of information technology to modern celestial 
navigation is very meaningful. This work is based on the 
concept of the genetic algorithm from the field of artificial 
intelligence. It combines a search method that mimics the 
natural process of biological evolution with the celestial fix 
principle of circles of equal altitude to rapidly and directly 
obtain accurate optimum vessel position data. Test results 
show that in addition to avoiding tedious and complicated 
computation and graphical procedures, this method is also 
more flexible, and is consistent with the development trend of 
celestial navigation in the electronic and information age. 


DEVELOPMENT OF CONTEMPORARY 
CELESTIAL NAVIGATION TECHNOLOGY 


The objective of celestial navigation, as traditionally 
practiced, is to determine the latitude and longitude of a vessel 
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ata specific time using observations of the altitudes of celestial 
bodies. This objective and the basic principles still remain 
unchanged today, and have formed a foundation on which 
a number of methods have been developed. The following 
is areview of the development of contemporary celestial 
navigation technology. 


Intercept Method 


The intercept method, which is the oldest method, was 
developed in 1875 by Marcq de St. Hilaire. This method 
translates each celestial altitude observation into a line of 
position (LOP) on the surface of the Earth. In principle, a series 
of observations defines a group of intersecting LOPs and this 
intersection represents the observer’s position the fix. Having 
been practiced for more than a hundred years, the intercept 
method is able to provide complete solutions for problems of 
astronomical vessel positioning. Its applicability continues 
in today’s maritime education and training, and in merchant 
shipping practice. However, it is an approximate method and 
involves graphical procedures. The entire process consists of 
observing, calculating, and plotting. Even when performed by 
professional seafarers, one astronomical positioning will take 
around 20 minutes. Considering the high speed of modern 
vessels, this method fails to achieve real-time positioning, and 
the lag time in positioning cannot guarantee navigation safety. 
Moreover, there are two theoretical assumptions implied in 
the intercept method. One is that the azimuth at the observer’s 
position with respect to the celestial body’s geographical 
position is the same as the azimuth at the assumed position. 
The other assumption is that when the co-altitude is sufficiently 
large, the circle of position can be taken as line of position. In 
view of these assumptions, the accuracy of astronomical vessel 
positioning is subject to the following two restrictions: 

1) As the choice of assumed position (AP) is by a trial-and- 
error method, its distance to the actual vessel position 
largely affects the accuracy of the result. Therefore, this 
distance should not exceed 30 NM. 

The altitude of the celestial body should not exceed 70 
degrees. Otherwise, the resulting error of curvature will 
increase when using line of position in lieu of circle of 
position on the Mercator Chart due to the reduced co- 
altitude. 


2 


Na 


Following the widespread application of information 
technology, some researchers applied the intercept method 
using computer programs. In Van Allen [13], the development 
is based on the geometry of the two LOPs. A least-squares 
approach for a multi-star fix is presented in Dewit [4] based 
on the plane geometry and straight lines formed by LOPs near 
the estimated position. This method is a direct mathematical 
translation of chart-based navigation. Although the intercept 
method has been computerized, the inherent drawbacks due 
to the two aforementioned assumptions still exist and limit 
its accuracy. 


Spherical Triangle Method (STM) 
and Vector-Matrix Method 


With the aid of computer programs, problems that could 
not be solved previously using the inspection table can now 
be solved using the STM or Vector-Matrix methods. With 
these tools, the celestial fix becomes fast and accurate. Unlike 
with the intercept method, there is no limitation on accuracy. 
Either the assumed position is not needed or the knowledge of 
an approximate position is sufficient. Some researchers have 
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studied a two-body fix [2, 3, 5, 6, 12, 14]. It is proposed that at 
a given time, the astronomical vessel position can be determined 
by observing circles of equal altitude of two celestial bodies. In 
a multi-body celestial fix, a statistical method such as the least 
mean squares technique is applied on the results of a number 
of two-body fixes to deal with the over-determined celestial 
fix problem and to solve for the final position [1]. However, in 
practice, it is rare to use two celestial bodies to obtain a vessel 
fix, and the statistical computation adds complexity to this 
method. Some other researchers have studied the multi-body 
celestial fix by invoking the circle of equal altitude. In this 
method, a celestial body’s altitude is viewed as a function 
of latitude and longitude, and a multi-body celestial fix is 
carried out directly, using the vector-matrix method, instead of 
indirectly using results of a two-body fix [9, 11, 15]. The vector- 
matrix formulations, while elegant, have the disadvantage that 
the solution still assumes prior knowledge of the observer’s 
position. Furthermore, the algorithm that proceeds from some 
initial value by sliding down the gradient of the goodness-of- 
fit parameter can converge to a local minimum that is not the 
best solution [7]. In addition, a multi-body celestial fix that is 
based on the least squares fit requires more than three celestial 
bodies. When there are only two celestial bodies available, 
this method cannot be used. In these cases, a virtual celestial 
body is introduced according to the vessel’s dead reckoning 
position. Consequently, the accuracy is affected by the assumed 
vessel position. 


GENETIC ALGORITHM TO SOLVE FOR 
ASTRONOMICAL VESSEL POSITION 


Advances in information technology have brought new 
changes to celestial navigation, especially in the early 1990s 
when research in this field blossomed. At this time, artificial 
intelligence began to be widely used in various research fields. 
Unfortunately, research on celestial navigation was soon 
largely replaced by the rapid emergence and development 
of GPS. At present, celestial navigation still serves as an 
important backup measure of vessel positioning in the ocean. 
It is therefore necessary to investigate the applications of 
information technology and artificial intelligence technology in 
celestial navigation and to improve the shortcomings of existing 
methods, which were described in the previous section. 


General principle of GA 


In recent years, the development of artificial intelligence 
technology has enjoyed mutual influence with the life sciences. 
The flourishing of evolutionary computation methods, such 
as the genetic algorithm (GA), reflects this characteristic and 
development trend. GAs are among the most effective methods 
for optimization problems. All types of objective functions 
and constraints can be processed without many mathematical 
prerequisites. During the search process, knowledge of the 
search space can be automatically obtained and accumulated. 
Accordingly, the search is modified to solve for the optimal 
solution or a less suboptimal solution. This method is very 
efficient and enables a parallel search. Two main features of 
GAs are a search strategy in a population and the exchange of 
information between individuals in a population. It simulates the 
overall learning process involved in the creation of a population 
from individuals. Starting from an initial population, in which 
each individual is a candidate solution in the search space of an 
optimization problem, the fitness of each individual is evaluated 
based on a fitness function. Multiple individuals with better 
fitness are subsequently selected, and then modified through 


genetic manipulations such as crossover (analogous to the 
exchange of information between individuals in nature) and 
mutation (where an individual is subjected to a sudden change 
and escapes the local optimum) to form a new population, thus 
achieving screening and evolution. The process is repeatedly 
executed. During the iteration process, every new population 
is superior to the population in the previous generation. 
Individuals continuously evolve toward better solutions, and 
eventually reach the optimum solution. 

There is always complexity and uncertainty associated 
with a celestial fix, which makes it difficult to obtain the exact 
solution mathematically. In this study, the superior search 
ability of GA is combined with fundamental positioning 
principles in celestial navigation to perform a fitting calculation 
for the optimum vessel position. 


The celestial fix principle of circles of equal 
altitude 


The central idea behind the celestial circle-of-equal-altitude 
fix is to find the best fit to the altitude of a celestial body 
observed as a function of time. The calculation altitude (Hc) 
of a celestial body is given as a function of the declination of 
body (Dec), the Greenwich hour angle (GHA), the observer’s 
assumed longitude (A), and the observer’s assumed latitude 
(L) by: 

Sin He; = SinL x SinDec;+ f 
+ CosLx CosDec; x Co(GHA; — À) a 


where i= | to n represents the celestial body and the observation 
data of the i celestial body 


The GHA and the Dec can be found from the time of 
observation using the Nautical Almanac or the electronic 
edition of Nautical Almanac. Hc in Equation | is the calculated 
altitude, which differs from the actual observed altitude. The 
altitude of a celestial body (hs) is measured by a sextant. An 
instrumental error correction (I), index correction (IC), and 
correction for optical deviations such as dip and refractions are 
applied to the sextant altitude. The resulting value is the actual 
observed altitude (Ho). Since the actual vessel position can be 
seen as a function of the altitude of the celestial body, which 
is a nonlinear function of L and A, we cannot solve equation 
1 directly. Therefore, the latitude and longitude of the actual 
vessel position may be obtained from the difference between 
Hc and Ho using an appropriate fitting algorithm. With the 
appropriate Dec and GHA in equation 1, any combination of 
altitudes of one or several celestial bodies can be used, and no 
plots or tables are required. Furthermore, when performing 
the multi-body celestial fix, different observing times lead to 
inconsistency in the observed altitude. It is therefore necessary 
to make corrections on the running fix. In this study, the 
correction approach proposed by Kaplan [8] is adopted. 


Implementation of a GA for the celestial fix 
problem 


The optimum vessel position in a majority of multi-body 
celestial fix problems is found by data fitting using the least 
mean squares method, and it is classified as an optimization 
problem. However, there still exists the problem of converging 
to a local minimum. GA shows excellent performance in 
optimum fix processing problems. Its main advantages include: 
(a) simple calculation without any complicated mathematical 
procedures; (b) it does not suffer from possible sensitivity to the 
initial position chosen to start the iterative process, meaning it 


has less of a tendency to converge toward the local optimum; 
(c) more flexibility in practice and rapid convergence toward 
the final result. The GA approach can be used to obtain a fix 
from multiple celestial observations of a single celestial body 
or of multiple celestial bodies (more observation data points 
produce better results). This method is in line with the needs 
of modern navigation. The procedure includes: 

1) Design of the objective function: The first step in solving an 
optimization problem is the design ofthe objective function 
to evaluate the fitness of each individual. A GA is used in 
this study to minimize the root mean squares of the altitude 
residuals (equation 2, root mean square error — RMSE), 
defined as the difference between each observed altitude 
(Ho) and the altitude computed from equation 1 (Hc). 
The residual defined differs from the altitude difference 
normally used to plot LOPs since the altitude difference 
is referred to an assumed position chosen to permit table 
lookup of whole degrees of Local Hour Angle. When 
the RMSE converges to meet the acceptance criteria, the 
optimum vessel position is obtained. 


(2) 


2) Design and encoding of individuals and population: In 
order to improve the efficiency and accuracy in this study, 
individual encodings in the GA are not in the form of 
traditional binary strings. Instead, real number encodings 
are adopted. Each individual contains two real number 
columns, representing latitude and longitude, to represent 
a possible solution (vessel position). Each population 
consists of multiple individuals. In generating the initial 
population, a reference position is introduced to improve 
search efficiency. This reference position is different from 
the assumed position in the intercept method. It is an 
approximate position of reference, which is not required 
to be within 30NM of the actual vessel position; indeed, 
it can be hundreds of nautical miles away. This permits 
a large degree of flexibility in choosing the reference 
position. This position is mainly a reference for a heuristic 
search. However, it is certain that a reference position closer 
to the actual position results in a smaller search space, 
resulting in higher efficiency and a better quality solution. 
The introduction of a new reference position does not add 
much inconvenience to seafarers. As stated by Pepperday 
[10], “If you know your approximate position. Why not 
use it?”. In fact, methods that claim that no DR reference 
position is required still need to assume prior knowledge 
of the observer’s position, or that it be obtained through 
complicated computational procedures. Moreover, DR 
position is also required for applying corrections on the 
running fix. The reference position is used as the datum 
point in this work. In the user defined maximum search 
degree and numbers of generations, random numbers are 
used to generate the initial position data of each individual 
in the initial population as the starting point of the evolution 
(equation 3, 4, 5): 


ae _ {double latitude 
Data structure of each individual Di: 3 
double longitude 


ee es IniL+Tol x RandNum 
Data of every initial individual (4) 
Ini +Tol x RandNum 
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where IniL and Inia represent the latitude and longitude values 
of the reference position, respectively. Tol is the user defined 
maximum search degree. RandNum is random real number 
between -1 and +1, and the random positions are generated to 
correct for the reference position. P, is the initial population, 
and n is the number of individuals in a population. 

3) Selection: The traditional roulette method is adopted as the 
selection method. This is a proportional selection strategy. 
In this method, individuals with a higher fitness level 
have a better chance of being selected to generate the new 
population. This concept is used in GA to mimic biological 
evolution. 

4) Crossover: The crossover rate is denoted as P.. n individuals 
in the population are stochastically paired to form n/2 
pairs. In order to determine in which pairs crossover takes 
place, the following procedure is repeated from i = 1 to 
n/2: Generate a random number r in the interval [0, 1], and 
if r<P., then apply crossover on the i pair. The crossover 
proceeds according to equation 6 and 7: 


{Lew =L1 +|L, -L2|x RandNum 


Diegi 
NOW" A New =21 +|; — Az |x RandNum 


(6) 


or: 
{Lew =L2 +|L, — L>|x RandNum 
Ane = Aq +A; — A, |x RandNum 


where D,,,, represents new individuals generated through 
crossover. L,,,, and Ay, denotes the latitude and longitude 
of new individuals generated through cross. L,, à, L,, A, 
are latitude and longitude of the parent individuals on 
which crossover is applied. RandNum is random real 
number between -1 and +1, and they are used to generate 
more optimum random positions. The crossover step 
can bring new positions closer (as shown in Figure 1) 
or further away (as shown in Figure 2). Since the parent 
individuals are already individuals with better fitness, it 
is more likely for the newly generated individuals to have 
further improvement (more closer to the actual vessel 
position). 


Dy (7) 


<)>Parent vessel position 1 


Vessel position after crossover 


€ Actual position 


Parent vessel position 2 


Fig. 1. Position becomes closer after crossover 


na Actual vessel position 


* Vessel position after crossover 


Parent vessel position 1 
Parent vessel position 2 


Fig. 2. Position becomes further away after crossover 
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5) Mutation: The mutation rate is defined as P,, The following 
procedure is repeated on individual i = | to n: Generate 
a random number r in the interval [0, 1], and ifr<P,, then 
apply the mutation operation according to: 


E Loig + Noise x RandNum 
or (8) 


A Mutation = “org + Noise x RandNum 


Mutation — 


D 


Mutation * 


D represents the data of individuals after mutation. 
L Mutation 200 Anuario, represent the latitude and longitude 
of individuals after mutation. Lo and Ao), represent the 
latitude and longitude of individuals before the mutation. 
Noise is a user defined degree of allowable perturbation, in 
place for starting new research. RandNum is random real 
number between -1 and +1, and these numbers are used to 
generate the stochastic perturbation degree. The mutation 
operation can prevent the search from converging to the 
local optimum (as shown in Figure 3). 
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Fig. 3. Mutation prevents the search from converging to the local optimum 


6) Improvement of the efficiency of the GA: In some cases, 
individuals in the population become increasingly similar to 
each other during the evolutionary process, thus limiting the 
GA’s global search ability, and causing the GA to converge 
toward a local optimum instead of the global optimum. To 
avoid this, we have enhanced the global search ability of 
the GA in the current study. At the early stage of evolution, 
a large crossover rate and mutation rate are imposed in 
order to maintain diversity in the population. At the later 
stages of evolution, individuals have already developed to 
a near optimal state. At this time, a large crossover rate or 
large mutation rate are likely to impair individual fitness. 
Therefore, the crossover rate and mutation rate are reduced 
at the later stages. 


Procedures of implementation 


Although a large population size is helpful in the search of 
vessel position, it inevitably increases the search time, making 
1t difficult to use for real time positioning. In the current study, 
the population size is set at 50 individuals. The crossover 
rate and mutation rate are 0.6 and 0.05, respectively. The 
maximum search space is set to 2 degrees (1.e. 4 degrees in 
latitude and longitude), which is sufficient for navigation. The 
perturbation range is set to 1 degree. Applying these settings 
in the genetic operations creates new generations to replace 
the parent generation. The search is terminated when either the 
RMSE converges to the acceptance criteria, or a fixed number 
of generations is reached. Termination of the search indicates 
that the optimal vessel position has been obtained. Experiments 
show that convergence to the accurate position generally occurs 
within 50 generations, and within approximately 5 seconds. 


generation = Max ? 


selection, crossover, mutation 
Selection 
Obtain new generation 


generation = generation + 1 


Fig. 4. The flowchart of GA for obtaining the optimum vessel position 


This satisfies the requirement of real-time positioning. The 


step-by-step procedures are (refer to Figure 4): 


1) Use real numbers to encode individuals in the solution 
space. Each individual corresponds to a candidate vessel 


position. 


2) Set the initial population size to contain 50 individuals, 1.e. 


stochastically generate candidate vessel 


3) 


4) 
5) 
6) 

Choose optimum 
vessel position 7) 
E 8) 
9) 


Calculate the value of the fitness function of each individual 
according to equation (1) and (2), and determine which 
individuals in the population are relatively close to the 
actual vessel position. 

Select individuals close to the actual vessel positions using 
the roulette technique. 

Apply crossover to the population. Crossover between 
vessel positions with good fitness generates vessel positions 
that are closer to the optimal positions. 

Apply genetic mutation to the population to avoid converging 
to local optima instead of the global optimum. 
Re-calculate the value of the fitness function of each 
individual according to equations (1) and (2). 

Determine whether the termination criteria has been met. 
If yes, then the genetic operation ends and proceeds to step 
(9), where the current optimum vessel position is exported. 
Otherwise, move back to step (4). 

Export the optimum vessel position. 


RESULT VALIDATION 


The validity of this study is verified by comparing results 


of a two celestial body fix and a multiple celestial body fix 
with those obtained by other methods. Visual Basic.Net 2010 


is used as the development tool. In order to add visual effects, 


components of Geographic Information System (GIS) are also 
introduced. This is also aimed at future celestial navigation 


education and integration with an Electronic Chart Display 


stochastically generate 50 candidate vessel positions. 
Tab. 1. Comparison of GA and relevant information from Hsu et al. [6] for two-body fix 


Body ZT 


Ho 


and Information System (ECDIS). 


GHA Dec 


Capella 20-03-58 


15*19.3” 


131924.8” 4558.4” N 


Alkaid 20-02-56 


77°34.9° 


003°14.2’ 49°25.7 N 


DR (20-04-00) 
Intercept Method 


L=41°34.8 N  1=017%00.5” W 
L=41°38.6 N  1=017%08.1? W 


L=41°39.17N à= 017°07.3 W 


Dec(d) Decim) Dech) DeeNS) GHA(d) GHA(m) GHA() Hr Mn See Hog) Holm) Hot) 


ss 


248 2 0 5% 15 193 | 
n Mus n n9 


17,1253504,41 654099, 0 002546 
171261,41 Seale O OEN 


76, 
17.122671,41 654176, 0.001450 


| GA Computston | 


Fig. 5. GA of two-body celestial fix 


POLISH MARITIME RESEARCH, No 3/2012 57 


Case 1 (2 Stars) 
--- Case 2 (4 Stars) 


1 11 21 31 41 
Generation 
Fig. 6. The evolutionary process for calculation of a fix 


Case 1 - Two-Body Fix: The data in this case study is taken 
from Hsu et al. [6]. The authors used the Simultaneous Equal- 
altitude Equation Method (SEEM) to perform a two-body 
celestial fix. Table 1 contains relevant experiment data. 

It can be seen from Table 1 that the observed altitude of 
Alkaid is as high as 77°34.9’, which exceeds the upper limit of 
70° of the intercept method. Hsu et al. [6] found that the LOP 
drawn using the intercept method shifted due to curvature error, 
thus resulting in an inaccurate vessel position estimate. 


The present study however is free from this drawback. In 
our experiments, even if the distance between the reference 
position and the actual vessel position exceeds 30NM and the 
search degree reaches 4°x4° (refer to Figure 5), the correct 
vessel position can still be obtained. Moreover, the RMSE 
converges to 0.0001 within 30 generations and within 3 seconds 
(refer to Case 1 in Figure 6). The calculation is performed in 
a rapid and accurate manner without the need to construct 
a visual celestial body. 

Case 2 - Multi-body Fix: There are four celestial bodies 
involved in this case study and more realistic in real marine 
navigation situation. Table 2 presents relevant observational 
data. In addition to the increased number of celestial bodies, 
correction on the running fix is also applied. The parameters 
in this experiment are in general identical to those in the 
previous case. However, the reference position is deliberately 
placed further so that it is almost 2 degrees away from the 
actual vessel position. Though the final RMSE (about 0.002) is 
slightly bigger than Case 1, caused by observation errors from 
multi-star sights, it can converge within 50 generations (refer 
to Case 2 in Figure 6). Results of this experiment show a good 
agreement with those obtained by the vector-matrix method or 
the graphical intercept method (refer to Figure 7 and Figure 
8). No plots or sight reduction tables are needed. The running 
time is within 5 seconds, and can meet the demands of real-time 
navigation fix. The applicability of the GA approach for both 
the two-body fix and multi-body fix is thus verified. 


Tab. 2. Multi-body fix 


Body ZT (1993/9/13) 


Ho 


GHA Dec 


Altair 18-00-00 


37°53.0° 


325%06.6” 08°51.£ N 


Fomalhaut 18-04-00 


27°54.0° 


279°24.2’ 29°39.1°S 


Achernar 18-08-00 


17°46.5’ 


240°21.7° 57°15.8 S 


Rasalhague 18-12-00 


41°35.5? 


002°04.8’ 12°34.1°N 


DR 


L=35°13.0 S 21=005%20.0”E Course = 220° Speed = 18 kts 


Intercept Method 


L=35°19.0°S 1=005%26.5” E 


Vector-Matrix 


2 Genetic Algorithm for Multi-Star Fix Astronomi 
Center of Reference Amea 


L=35"18.6"S 1=005%26.9 E 
L=35"18.6"S  1=005%27.0” E 


66 


27 


.428925,35 315829, 0.010731 
28325,35 315829, 0.010731 
}28825,35 315829, 0.010731 


| GA Computation 


Fig. 7. GA of multi-body celestial fix 
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Fig. 8. Comparison of GA and intercept method in multi-body celestial fix 


CONCLUSIONS 


Although celestial navigation has its shortcomings, it cannot 
be replaced by GPS due to the advantages of low cost, high 
reliability, strong independence, and inability to be detected. In 
this study, a genetic algorithm that mimics biological evolution, 
taken from the field of artificial intelligence, is applied to the 
celestial navigation fix problem. The main features of the 
GA include: (1) simple calculation without any complicated 
mathematical procedures; (2) Avoidance of converging toward 
a local optima; (3) more flexibility in practice and rapid 
convergence toward the final position estimate. The algorithm 
can be used to obtain a fix from observation of a single celestial 
body (multiple observations) or of multiple celestial bodies, and 
more observation data produces better results. Experimental 
results show that in addition to evading tedious and complicated 
computation and graphical procedures, the GA approach is 
more flexible compared with other methods. This approach also 
has the potential to be integrated with ECDIS. Moreover, it can 
also be used to confirm whether GPS is functioning normally. 
Navigation technology is moving into the age of e-Navigation, 
and the continuous development of celestial navigation can only 
be realized through integration with information technology. 
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Automatic Identification System (AIS) 
application to anti-collision manoeuvring 
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ABSTRACT 


Common use of shipboard AIS creates conditions for the use of a new kind of dynamic data in the situation 

of the risk of collision. AIS position report is a source of supplementary information derived from error 

leveraged radar measurement. However, in view of the results of the studies there are opinions with regard 

to inconsistent AIS dynamic data in the process of decision-making by the officer of the watch. By taking 

into consideration the recordings of the studies and technical specification of AIS it can be concluded that 
the results of inconsistent data have significant role in collision avoidance manoeuvring. 


Key words: AIS; collision avoidance 


INTRODUCTION 


According to International Regulations for Preventing 
Collisions at Sea, Rule V “Look-out” [7] „Every vessel shall 
at times maintain a proper look-out by sight and hearing as 
well as by all available means appropriate in the prevailing 
circumstances and conditions so as to make a full appraisal 
of the situation and of the risk of collision”. Effective look-out 
should ensure early detection of objects. Another problem is to 
assess ship motion parameters and foresight situation. One of 
the effective undertakings is radar look-out in any conditions 
of visibility, especially, when the vessel navigates in area 
of intensity traffic or at night. It is common knowledge that 
Radar and ARPA have some efficacy limitations [8]. Detection 
of small objects is limited by sea clutters and unfavourable 
weather conditions (rain, snowstorm), and by close frequency 
working radio-transmitters. 

Characteristic phenomenon is death-zone appearance 
(minimum range of detection), radar shadow effect (there are 
reflections from funnel, mast and other constructions on the 
vessel) and wave reflection. The radar information, usually 
presented in polar coordinates system is subject to some 
measurement error exceeding approximately 1% of radar range. 
As far as ARPA devices are concerned, automatic tracking 
objects can be lost during torrential ship manoeuvring or when 
passing ships within close distance. Therefore it is reasonable to 
supplement the information derived from radar with additional 
information delivered automatically and continuously with 
better accuracy. A system which appears suitable to provide 
this kind of information is AIS which may be classified as radio 
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navigation system which uses radio waves to transmit data with 
regard to the ship motion parameters. This information is more 
accurate than radar information, however for many sea officers 
it looks less reliably than own radar information. 


QUALITY OF AIS INFORMATION 


To illustrate AIS information with respect to the movement 
of objects in the operation area of the system the following 
elements of the ,,Position Report” containing the AIS messages 
1, 2 and 3, namely: 

- SOG,(Speed over ground); 
- ROT,(Rate of turn); 
- HDG, (Heading); 
- geographic position; 
- time stamp; 
can be used. 


Dynamic information has been used in the calculations 
performed by MADSS (Multi-agents Decision Support 
System), described by [2]. On the basic of received messages 
AIS calculates ship motion parameters and generates new 
motion parameters of own vessel (course, speed), which leads 
to pass by one another depending on calculated CPA (Closest 
Point of Approach). 

Common use of AIS caused the emergence of the opinion 
on the imperfections of the system related to the lack of 
transmission or the transfer of not reliable information. 
Studies on the incompleteness and integrity of AIS information 
published to date are generally linked to the message No. 5 
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20% 
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16% 
14% 
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NAV.STAT SOG 
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[HEAD 


COG ROT 


N 


IMSI C.SIGN S.TYPE 
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A % of Incompleted messages 
EE % of vessels with incompleted messages 
Fig. 1. Initial results of incomplete AIS information [acc. these authors’ studies]. Abbreviations: NAV.STAT — Navigational status; POS — Coordinates of 
ship s position; SOG — Speed over ground; COG — Course over ground; T.HEAD - true heading; ROT — Rate of turn; MMSI - Maritime Mobile Service 
Identity number of the ship; IMO — Ship 5 IMO number; C.SIGN — Ship 5 call signal; S.NAME — Ship s name; S. TYPE — Type of the ship; 
M.DRA — Mean draft of the ship; S.DEST — Destination of the ship; 


(Static & voyage related data). This problem was also analysed 
by [4] and [1, 3, 5, 6]. However, it seems that there are no 
publications concerning quality of the dynamic information, 
particularly the one which is important in the analysis of anti- 
collision maneuvering. In addition mentioned publications has 
been prepared in the beginning stage of AIS. In view of this 
facts, an initial investigation has been conducted by authors to 
determine to what extent the information transmitted by AIS 
and derived from ship sensors is complete. 


INITIAL INVESTIGATION OF THE 
AIS DYNAMIC INFORMATION 
INCOMPLETENESS 


Completeness of the messages No. 1, 2, 3 based on analysis 
of contents of the messages AIVDM has been studied. Analysed 
data originate from the Gulf of Gdansk and contain recorded 
information of AIVDM mnemonic from 2006.04.12, at 00.00 
to 2006.04.12, at 23.59. The results of the studies are presented 
in Figure 1. Dark grey bar shows the percentage of messages 
which contain incomplete information, and the light grey bar 
shows the percentage of vessels responsible for this state. 

The AIS information studies were based on recorded 
messages in text files containing the messages (received by 
AIS) about vessels located in the Gulf of Gdansk during one day 
only. Preliminary analysis of the results leads to the conclusion 
that the biggest indicator of unfitness is characterised by 
information about the rate of turn which seems to be very 
important for anticollision activity. Incompletness of messages 
is at the level of 23%. These messages were sent by 24% of 
ships, approximately. Additionally, it is worth to noting the 
incompleteness of information concerning True Heading, which 


amounts about 6% and was observed in the case of 19% of 
vessels. Therefore, there is a reason to perform detailed studies 
on the incomplete information concerning True Heading and 
Rate of Turn. 


STUDIES OF INCOMPLETENESS OF 
SELECTED ELEMENTS OF AIS DYNAMIC 
INFORMATION 


On the basis of technical specification, AIS correct statement 
of data concerning volume and value of the incomplete 
information are shown in Tab. 1. 

According to the assumptions, detailed studies of dynamic 
information has been performed in 2010 and 2011 in selected 
49 days. Data has been recorded with shore receiver situated in 
the building of Polish Naval Academy in Gdynia. Information 
about True Heading and Rate of Turn characterise the biggest 
information unfitness like in 2006 preliminary study. Summary 
of this are presented in Tab. 2. 


Tab. 1. Summary of the ranges of correct 
and incorrect data in the message No 1 


Correctly value 


Type of 
information 


Incomplete 
information 


511 
(149 hex) 


128 
(80 hex) 


TRUE HEADING 


[0.359] 


RATE OF TURN 


127 


Figure 2 presents summary of the studies concerning 
AIS information unfitness, and takes into consideration 
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Fig. 2. Summary of unfitness studies of AIS information concerning Rate of turn [acc. 


these authors’ studies] 
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Tab. 2. Summary of unfitness studies concerning True heading and Rate of turn [acc. these authors’ studies] 


TRUE HEADING [%] 


RATE OF TURN [%] 


incompleteness 
(AIVDM sentences) 


2006.04.03 


unfitness 
(ships) 


incompleteness 
(AIVDM sentences) 


unfitness 
(ships) 


2006.04.04 


2006.04.05 


2006.04.06 


2006.04.07 


2006.04.08 


2006.04.10 


2006.04.11 


2006.04.12 


2006.04.13 


2006.04.18 


2006.04.19 


2006.04.20 


2006.04.21 


2006.11.22 


2006.11.23 
2006.11.24 


2006.11.25 
2006.11.26 


2006.11.27 


2006.11.28 


2006.11.29 


2010.12.16 


2010.12.17 


2010.12.26 


2010.12.27 


2010.12.28 


2010.12.29 


2010.12.30 


2010.12.31 


2011.01.01 


2011.01.02 


2011.01.03 


2011.01.04 


2011.01.05 


2011.01.06 


2011.01.07 


2011.01.08 


2011.01.09 


2011.01.10 


2011.01.11 


2011.01.12 


2011.01.14 


2011.01.15 


2011.01.16 


2011.01.17 


2011.01.18 


2011.01.19 


2011.01.20 


2011.01.21 


POLISH MARITIME RESEARCH, No 3/2012 


30% 


25% 


ANNAN NANNY 


20% |; 


AAN ANNAN NAS 


SANA, 


A 


343% 


ANANNENANNANNANA NAAA 
O 


ANNAN 
ANNAN 


A 
d AA 
A 


15% 


10% 


5% 


0% 


2006.11.28 eneee 
2010.12.16 Free 


2010.12.26 
2010.12.28 
2010.12.30 
2011.01.01 
2021.01.03 


2006.04.07 Rm 
2021.01.05 


2006.04.03 
2006.04.10 


2006.04.05 
2006.04.20 


2006.04.12 
2006.04.18 
2006.11.22 
2006.11.24 
2006.11.26 


SNA 


> N 
AAO 


2021.01.07 
2021.01.20 p 


2021.01.09 
2021.01.11 7 
2021.01.14 Je 
2021.01.16 
2021.01.18 


Fig. 3. Summary of the studies on unfitness of AIS information concerning True Heading [acc. these authors’ studies] 


incompleteness of Rate of Turn. The Dark grey bar presents 
percentage of incomplete messages concerning Rate of Turn, 
whereas the light grey bar presents percentage of ships, which 
are responsible for this state. 

Figure 3 presents results of the incompleteness studies 
of AIS information concerning True Heading. As previously 
mentioned, dark grey bar presents percentage of incomplete 
messages, whereas the light grey bar presents percentage of 
ships, which are responsible for this state. 

Taking into consideration the investigations performed 
in 2010 and 2011, only an imperceptible decrease of HDG 
and ROT coefficients with “AIVDM sentences” criterion is 
observed. The analysis of the results is presented in Tab. 3. 


Tab 3. Analysis of the studies on AIS information unfitness, gained in the 
period from 2006.04.03 to 2006.11.29 


RATE OF TURN 


Sentences 
(AIVDM) 


24.70% 
24.85% 
32.70% 
16.81% 
3.68% 
0.14% 


TRUE HEADING 


Sentences 
(AIVDM) 


18.27% 
17.89% 
22.37% 
14.61% 
2.25% 
0.05% 


ships 


Ships 


21.74% 
20.97% 
26.67% 
17.65% 
2.46% 
0.06% 


m, - median, x - arithmetic mean, 
o - standard deviation, o? - variance 


Tab. 4. Analysis of the studies on AIS information unfitness, gained in the 
period from 2010.12.16 to 2011.01.21 


RATE OF TURN 


Sentences 
(AIVDM) 


13.19% 
13.48% 
17.75% 
7.80% 
2.42% 


TRUE HEADING 


Sentences 
(AIVDM) 


13.58% 
14.42% 
17.38% 
7.69% 
2.59% 


ships Ships 


23.36% 
23.60% 
27.62% 
18.56% 
2.21% 


24.29% 
24.56% 
28.57% 
18.56% 
2.49% 


m, - median, x - arithmetic mean, 
o - standard deviation, o? - variance 


The essential percentage of misinformation about the 
Heading and Rate of Turn observed during investigations 
demands further explanations. One of reasons can be the 
frequent lack of measures of ROT on smaller ships. The 
other sources can be moored ships with active gyro, what 
can give the effect of misinformation about the Heading. 
Possible source of unfitted information about HDG can 
be some problems with connection between gyro and AIS 
device. Interesting is observation that different strategies of 
investigations gives us different results. For example mean 
value of wrong information about ROT is 13% or 23% - 
depends the name of ship is analyzed or type of AIS sentence. 
This leads to the conclusion that probably some ships without 
sources of ROT be present in analyzed radius more times 
than another. 


CONCLUSION 


- Potential application of AIS information for anti-collision 
manoeuvres seem to be promising however concerning 
True Heading and Rate of Turn data transmitted by this 
system rises a number of questions among practitioners. 
The question is why so many times information about 
HDG is unfitted. Suggestion is that there are some 
technical problems with connection between gyro and 
AIS device. 

- During presented studies two approaches have been 
applied. The unfitness rate was based on analysis of 
messages in the system and on the number of ships 
transmitting this messages. The arithmetic means of 
statistical features take similar values. The variance 
of data reveals their high variability characteristics of 
the study and the incompleteness of Rate of turn is at 
the level of 15.89% (Tab. 3). The variance of random 
variables provides information about their low dispersion. 
It is worth noting that the sample taken for testing was 
gained only from the Gulf of Gdansk, and the number 
of statistical units gives only a general view of the 
information unfitness during anti-collision manoeuvring. 
Conclusion about possible using AIS information for 
anti-collision tasks is that more reliable is information 
about course over ground transmitted by the ship than 
her heading. 
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ABSTRACT 


Over the last decades a really tough competitive environment has been observed in shipbuilding sector. Under 

these circumstances shipyards desire to reduce product cycle time in order to manufacture the product as 

soon as possible. For this reason it is required to make some alterations in shipyard production system in 

order to reduce cycle time. In this study a panel fabrication line used in a shipyard, has been considered. 

The workstations of the panel line have been modeled by using ARENA simulation software and the effects 
of matrix module assembly on panel line throughput, have been determined. 


Key words: shipyard; ship production; panel line; matrix module assembly; simulation 


INTRODUCTION 


In recent years there is a hard competitive environment in 
many industry branches including shipyard industry. In order 
to increase their competitiveness shipyards have to consider the 
factors affecting competitiveness according to Rashwan and 
Naguib [1]. One of the factors is delivery time. To be able to 
keep the competitive power, the companies have to deliver the 
products to their customers on time and reduce the production 
cycle time. For this reason the companies have to investigate 
their production systems and improve their processes. In 
order to see the effects of the improvements on a real system, 
simulation tool is mostly used. 

In this study simulation is used as an optimization tool. 
Simulation is used to understand character of real systems. There 
are many simulation applications in shipbuilding, as described 
in the literature. In a common study realized by Michigan 
University and Seoul National University, the whole processes 
ofa shipyard are attempted to be modeled with simulation and 
the effects of some changes on the system are perceived [2]. 
In the study of Okumoto ef al. [3], performed by modeling the 
scaffold placement with three-dimensional simulation CAD 
system, the effects of the changes on the system have been 
observed. Shin [4] has modelled the workstations of the sub- 
assembly line with the use of simulation and determined the 
effect of placing a welding robot on productivity. In the other 
study of Shin [5] an optimum shipbuilding layout has been 
found by using simulation. Alfeld et al. [6] has used a special 
software in simulating the shipyard processes to ensure that 
the planners take decisions easily. Alkaner [7] has simulated 
the processes of the profile cutting station and by making 
some changes in resources, the effects of these changes have 


been investigated. Doing some alterations on panel production 
station, Greenwood ef al. [8] has investigated the effects of 
the changes on the production system. Lee ef al. [9] has found 
the effects of the intelligent welding robot system on welding 
performance. Cha and Roh [10] have developed a simulation 
framework and applied it to block erection process. As can be 
concluded, there are many application fields of simulation in 
shipyard industry. 

In this study the processes of the panel line of a shipyard 
located in Turkey have been considered. Firstly, the detailed 
process analysis of the panel line has been performed. In 
this way the whole work activities of the panel line and their 
durations have been found. In the second stage of the study 
the simulation model of the panel line has been elaborated by 
using ARENA 11.0 software. The required data achieved from 
the process analysis have been put in the simulation model. 
Then the model has been run for 10-day period and number of 
products manufactured by the panel line has been determined. 
In the third stage of the study an alteration in section spot 
welding station has been made. And, this alteration has been 
put in the simulation model and the model has been run again 
for 10- day period. Then the effect of this alteration on the 
system throughput has been determined. 

As mentioned above, the alteration made in the panel line 
is interesting in the case of section spot welding station. In 
the section spot welding station of the panel line, minor and 
sub- assemblies are mounted on the flat panel assembly and 
their matrix structure is formed. In many shipyards in Turkey, 
minor and sub-assemblies are welded one by one on flat panel 
assembly in order to produce matrix structure. In this way it 
takes a longer time to manufacture the matrix structure and 
a bottleneck situation may occur on the panel line. Instead of 
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the mounting of minor and sub-assemblies one by one, they 
should be mounted on flat panel assembly as a module. In 
other words, the matrix structure can be mounted as a module. 
It is believed that mounting the matrix structure as a module 
reduces the production cycle time and increases the throughput 
of the panel line. 


METHODOLOGY 


In this study, the steps presented in Fig 1 are followed, 
respectively. At the beginning the workstations forming the 
panel line are determined and defined (Step 1). Then the 
product to be manufactured in the panel line and its sub- 
compenents are defined (Step 2). After the process analysis of 
the workstations is performed, the panel line is modelled by 
using simulation (Model 1), (Step 4) and the model is run for 
a specified duration time (Step 5). In Step 6, work flow in the 
section spot welding is changed and module mounting is carried 
out instead of mounting the minor and sub-assemblies one by 
one on the flat panel assembly (Stage F). In Step 7, the panel 
line is modelled by using simulation again and consequently 
Model 2 is achieved. In Model 2, the workstations’ work flows 
remain constant except for the section spot welding. In Step 8, 
Model 2 is run for the same time period as in the case of Model 
1 and throughput quantity of the panel line is achieved. At the 
last step (Step 9), Model 1 and Model 2 are compared mutually 
in terms of throughput quantity. 


Determination of panel line workstations (Step 1) 


Definitions of product to be manufactured in panel line 
(Step 2) 


Process analysis of panel line workstations (Step 3) 
Modeling panel line by using simulation (Model 1) (Step 


Running simulation for a specified duration time and 
determination of panel line throughput (Step 5) 


Process analysis of section spot welding with module 
mounting (Step 6) 


Modeling panel line by using simulation (Model 2) and 


determination of panel line throughput (Step 7) 


Comparison of Model 1 with Model 2 (Step 8) 


Fig. 1. Methodology of the study 


Determination of panel line workstations (Step 1) 
The panel line is a production cell where flat structures 


are fabricated. It consists of different types of workstations. 
There are 9 workstations on the panel line. Each workstation 
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has a function in the production process. Tab. 1 shows the 
workstations located on the panel line. 


Tab. 1. Workstations on the panel line 


No. of 


Workstation name 
workstation 


Edge cutting station 


Edge cleaning and sequencing 


Panel production 


Panel cutting 


Profile spot welding 


Profile tig welding 


Section spot welding 


Section tig welding 


Grinding 


The edge cutting operations of the plates are performed 
in the edge cutting station (11). The cutting operation is 
carried out by using plasma. The plates being cut are then 
sent to the edge cleaning and sequencing station (12) where 
the edge cleaning operation of the cut surfaces is performed. 
Grinding machines are used in this operation. The sequencing 
operation is also carried out in this station in order to sequence 
the plates which enter the panel production machine. The 
plates are then sent to the panel production station from the 
edge cutting and sequencing station. In the panel production 
station (13), the plates are mounted by using submerged arc 
welding. As an output of this station a panel is produced and 
then it goes to the panel cutting station where inside and 
outside cutting operations of the panel are performed. In the 
panel cutting station (14), marking operation is also done. 
The cut and marked panel is sent to the profile spot welding 
station (15). Here, the profiles are mounted on the flat panel 
and a flat panel assembly is produced as an output. a spot 
welding machine is used for this operation. The flat panel 
assembly is then sent to the profile tig welding station (16) in 
order to complete the welding operations. After completion 
of the welding operation of the flat panel assembly it goes to 
the section spot welding station (17) where minor assemblies 
and sub - assemblies are mounted on the flat panel assembly 
by using spot welding. As an output, a major sub-assembly is 
manufactured in this station. The major sub-assembly is sent 
to the section tig welding station (18) from the section spot 
welding station in order to complete the welding operation. 
Finally, the major sub-assembly arrives at the grinding station, 
the last workstation of the panel line (19). In this station the 
grinding operations of the welded places of the major sub- 
assembly are carried out. After completion of the grinding 
operations it leaves the panel line. Fig. 2 shows the material 
flow on the panel line. 


Definitions of products to be manufactured in 
panel line (Step 2) 


In the panel line interim products which have flat 
structure are produced. For a double bottom block, flat panel 
assembly and major sub-assembly are fabricated in the panel 
line. In this study, a major sub-assembly was considered as 
a product. 

In ship production some codes, each representing 
a production stage of blocks, are used. Such coding system is 
very useful to seperate and check the production stages orderly. 
Tab. 2 shows the production stages and their definition. 


Section tig 
welding 


Profile tig 
welding 


Section spot 
welding 


Profile spot 
welding 


TT 


Panel 
production 


Edge cleaning and 
sequencing 


Edge cutting 


Fig. 2. Work flow through the panel line 


The major sub-assembly (Stage G) is regarded as 
a throughput of the panel line in question and it includes 
structures of various types, such as single section parts (Stage 
A), single plate parts (Stage B), minor assemblies (Stage C), 
sub-assemblies (Stage D), flat plate assembly (Stage E), and 
flat panel assembly (Stage F), as defined in Tab. 2. 

The single section parts and single plate parts which have 
specified dimensions are described as a and B production 
stages, respectively. a single profile is defined as a production 
stage and a single plate is defined as B production stage. 

If one single section part and one single plate part are 
assembled together, the minor assembly (C) is manufactured. 
If two or more minor assemblies are fitted together the sub- 
assembly (D) is built. 

The flat plates constitute flat panel structures. If two or more 
flat plates are fitted together they form the flat plate assembly 
(E). If single section parts (A) are fitted onto the panel, the panel 
with profiles, called the flat plane assembly (F), is formed. 

Minor and sub-assemblies (C and D production stages) are 
fitted onto the flat panel assembly (F) to form finally the major 
sub-assembly (G). 

As above mentioned, the major sub-assembly includes 
various types of production stages. Fig. 3 shows the product 
breakdown structure of the major sub-assembly considered in 
this study. 


Tab. 2. Production stages and definitions 


Production Detiniions Structures representing 
of production > 
Stage production stage 
stages 
— 


Single section 


A part 


Single plate 
part 


Minor 


= assembly 


Sub assembly 


Flat plate 
assembly 


Flat panel 
assembly 


Major sub 
assembly 


In the major sub-assembly production process some 
mounting operations take place. In the first stage a set of single 
section parts of 7 in number are welded together and the flat 
plate assembly is formed. In the second stage the flat plate 
assembly and a set of single section parts of 18 in number are 
welded together and the flat panel assembly is produced. In the 
third stage, the flat plate assembly, a set of minor assemblies 
of 14 in number and a set of sub- assemblies of 9 in number 
are mounted together to form the major sub-assembly being 
a throughput of the panel line. 


Process analysis of panel line workstations 
(Step 3) 


So far, the panel line and the structure of major sub- 
assembly as an output are briefly discussed. In this section 
a detailed process analysis of the panel line production 
system has been performed. During the process analysis 
each of work stations in the panel line has been considered 
in detail. 

The main point of the process analysis is to determine 
the work activities. After determining work activities their 
operation times are calculated. Then, by considering the parallel 
and serial work activities, the completion times of work stations 
are determined. It is impossible to present here all the work 
activites in panel line. For this reason the process analyses of 
the profile spot welding (15) and profile tig welding (16) stations, 
are only exemplified in Tab. 3 and 4. 
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Major sub-assembly 
(Production stage G) 


(1) 


Flat Panel Assembly Minor Assembly 
(Production Stage F) (Production Stage C) 
(1) (14) 


Sub-assembly 
(Production Stage D) 


(9) 


Flat Plate Assembly 
(Production Stage E) 
(1) 


Single Section Part 
(Production Stage A) 
(18) 


Single Plate Part 
(Production Stage B) 
(7) 


Fig. 3. Product breakdown structure of major sub-assembly 


Tab. 3. The process analysis of the profile spot welding station (I5) 


No. of ae. eae Activity duration 
oe Activity description f j 
activity time [min] 


The operator walks to the crane 


ji 


The operator runs the crane 


The crane goes to profile stock area 


The operator assistants go to profile stock area 


The crane comes down the profile 


The crane holds the profile 


The crane lifts the profile 18.037 


The crane transports the profile from profile stock area to the porter system 8.473 


SLD] Nn] BR] Ww] bo 


O| 0 


The workers walk to the porter system 3.609 


= 
o 


The crane puts down the profile onto the porter system 12.274 


rn 
hb 


The crane leaves the profile operation area 4.428 


= 
N 


The workers settle the profile on the porter system 3.8 


jà 
(0S) 


The operator walks to the proter system 0.118 


Eh 
& 


The workers walk to the profile welding area 0.404 


= 
Nn 


The operator runs the porter system 0.166 


n 
oN 


The operator drives the porter system to the welding area 2.926 


jà 
N 


The operator walks to profile spot welding machine 0.042 


— 
00 


The operator cleans the welding torch 1.5 


The operator runs the profile spot welding machine 0.5 


pă 
\o 


N 
o 


The profile spot welding machine goes to the porter system 44.755 


N 
he 


The profile spot welding machine comes down the profiles 3.8 


The profile spot welding machine transports the profile from the porter 


N 
N 


46.486 
system to the flat plate assembly 


The profile spot welding machine comes down the profile onto the flat plate 


N 
Ww 


; : 111.394 
assembly and alignment is performed 


The profile spot welding is prepared for welding operation. 


The process of spot welding 


The conveyor system transports the flat plane assembly 


Total activity duration time 
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Tab. 4. The process analysis of the profile tig welding station (16) 


Activity description 


The conveyor transports the flat panel assembly to the tig welding station 


Activity duration 


time [min] 


The operator removes the slag from the welding torch 


The operator drives the tig welding machine to the starting point of welding 


11.577 


The operator takes down the welding torches on the welding area 


17.1 


The process of tig welding 


301.071 


The operator takes up the welding torches 


3.154 


The conveyor transports the flat panel assembly to the buffer area 


4.134 


Total activity duration time 


413.61 


Tab. 5. Completion times of the workstations of the panel line (Model 1) 


No. of station Station name 


Edge cutting station 


Number of work activities 


Station completion time 
[min] 


TRIA(111.5,144.9,200.8) 


Edge cleaning and sequencing 


TRIA(119.2,154.9,214.6) 


Panel production 


TRIA(368.2,478.6,49 1.5) 


Panel cutting 


TRIA(226.1,295.5,409.1) 


Profile spot welding 


TRIA(174.05,227.5,233.3) 


Profile tig welding 


TRIA(212.7,279.3,386.7) 


Section spot welding 


TRIA(501,656.3,911.8) 


Section tig welding 


TRIA(278,361.4,506) 


OlaolrAL NI NAR] WwW] mfr 


Grinding 


Tab. 3 illustrates the work activities of the profile spot 
welding station. As can be seen from Tab. 3, there are 298 
work activities to be perform in the profile spot welding station 
and its total activity time amounts to 380 min; whereas the 
station completion time is only 227.5 min because some of the 
activities are parallel. That is why the total activity time and 
the station completion time are not same. Tab. 4 represents the 
work activities of the profile tig welding station. In this case 
110 work activities are performed during the total activity time 
of 413.61 min; whereas the station completion time amounts 
to 279.3 min. 

In the same way the process analysis of other workstations 
of the panel line are performed and the station completion times 
are achieved, as shown in Tab. 5. 

It should be noted that the station completion times 
calculated from a comprehensive process analysis are regarded 
as optimistic ones. Because their distribution is assumed 
triangular, expected and pessimistic times are also needed to 
be assigned. In this study the optimistic and pessimistic times 
are assigned from gained experience. 


Modeling panel line by using simulation 
(Model 1) (Step 4) 


In this step the ARENA 11.0 software has been used for 
modeling the panel line. The required data have been achieved 
from the process analysis (Step 3); on this basis, apart from 
station completion times, transportation duration times are 
achived. The duration times are calculated by considering 
the production system and they are thought to have triangular 
distribution. Tab. 6 shows the duration times of transportation 
between workstations. 


TRIA(85,111.3,154.7) 


Tab. 6. Duration times of transportation between workstations 


Between workstations | Transportation times [min | 


Arrival to I1 TRIA(1.7,2.2,3) 

n>1 TRIA(1.4,1.8,2.5) 
2-13 TRIA(2.6,3.3,4.6) 
13—14 TRIA(0.7,0.9,1.2) 
1415 TRIA(0.5,0.6,0.9) 
15—16 TRIA(1.2,1.5,2.1) 
16-17 TRIA(2.1,2.7,3.7) 
1718 TRIA(1.2,1.5,2.1) 
I8—19 TRIA(1.1,1.4,1.9) 

In the simulation model in question, machine failures are 

also taken into considerations to reflect the real environment, as 


shown in Tab. 7. These values are not calculated but estimated. 
Failure times are thought to have exponential distribution. 


Tab. 7. Failure times of workstations 


Down time 
[min] 
EXPO(10) 
EXPO(4) 
EXPO(15) 
EXPO(20) 
EXPO(12) 
EXPO(18) 
EXPO(5) 
EXPO(25) 
EXPO(8) 


Station name 


EXPO(120) 
Edge cleaning and sequencing | EXPO(140) 

Panel production EXPO(180) 
EXPO(160) 
EXPO(130) 
EXPO(155) 
EXPO(165) 
EXPO(200) 
EXPO(160) 


Edge cutting station 


Panel cutting 
Profile spot welding 


Profile tig welding 
Section spot welding 
Section tig welding 
Grinding 
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In Figure 4, the simulation model of the panel line in Tab. 8. Module definitions 


question is presented. Module 


Table 8 shows the module definitions of simulation model Module name ae Module name 


Part Arrival I5 Arrival Station 
Arrival Station Profile Spot Welding 
Route to I1 Station Route to I6 Station 
Il Arrival Station 16 Arrival Station 
Edge Cutting Profile Tig Welding 
Route to I2 Station Route to I7 Station 
12 Arrival Station 17 Arrival Station 
Edge Cleaning Section Spot Welding 
Route to I3 Station Route to 18 Station 
13 Arrival Station I8 Arrival Station 
Panel Production Section Tig Welding 
Route to 14 Station Route to I9 Station 
14 Arrival Station 19 Arrival Station 
Panel Cutting Grinding 
Route to I5 Station End of Panel Line 


in Fig. 4. 


DE Se EH 
bo Lor] fie Sar a 
PEE AHA eee 


BH Ha] 


Fig. 4. Simulation model of the panel line 


1 
2 
3 
4 
5 
6 
7 
8 
9 
10 
11 


Rh 
N 


= 
UY 


do 
= 


Tab. 9. The process analysis of the section spot welding station (17) (Model 2) 


Activity duration 
time [min] 


E Activity description 
Transportation of matrix module structure to flat panel assembly 
Alignment of matrix module structure on flat panel assembly 
Horizontal spot welding and grinding operation after spot welding 
Workers go to pick up single plate parts 
Workers put down single plate parts 
Workers do alignment of single plate parts 
Horizontal spot welding of single plate parts 
Vertical fixing of single plate parts 
Vertical spot welding of single plate parts 
Operating crane 
Crane goes to single plate parts stock area 
Crane comes down onto single plate part surface area. 
Crane holds single plate parts 
Crane lifts single plate parts 
Crane transports single plate parts to flat panel assembly 
Crane puts down the single plate parts for marking 
Horizontal fixing of single plate parts 
Crane leaves the surface area of single plate parts 
Horizontal spot welding of single plate parts 
Vertical fixing of single plate parts 
Vertical spot welding of single plate parts 
Operating grinding machine 
Vertical and horizontal grinding after spot welding 
Crane goes to pick up lifting lug 
Crane comes down onto lifting lug’s surface area 
Crane holds lifting lug 
Crane picks up lifting lug 
Crane transports lifting lug to flat panel assembly 
Crane puts down lifting lug on flat panel assembly 
Fixing lifting lug on panel 
Crane leaves lifting lug’s surface area 
Spot welding of lifting lugs 
Crane departs from lifting lug 
Cleaning 
Transportation of major sub-assembly 
Total duration time 
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Running simulation for a specified duration 
time and determination of panel line 
throughput (Step 5) 


The simulation model has been run for 10-day period under 
the assumption that the shipyard operates in two shifts. Each 
shift lasts 8 hours. Number of replication of the model is equal 
to 5. As a result of the running, the panel line has produced 11 
major sub-assemblies. 


Process analysis of section spot welding with 
application of module mounting (Step 6) 


When the matrix module structure is assembled on flat 
panel assembly, the completion time of the section spot welding 
station will obviously change by nature of the things. To see 
the effect of the changing on throughput, the new completion 
time of the section spot welding station should be put in the 
simulation model shown in Fig. 4. In the current panel line 
system in question the completion time of the section spot 
welding station is determined as shown in Tab. 5. Tab. 9 
presents the work flow of the matrix module structure. When 
the matrix module structure is assembled on the flat panel 
assembly the completion time of the section spot welding 
station reaches 284 min. 


Modeling panel line by using simulation 
(Model 2) and determination of panel line 
throughput (Step 7) 


In this step the simulation model shown in Fig. 4 is 
applied. All the completion duration times of the workstations 
remain constant except for the section spot welding station. 
By changing the completion time of the section spot welding 
station, Model 2 was obtained. The completion times of the 
work stations for Model 2 are shown in Tab. 10. 


Tab. 10. Completion times of the workstations on the panel line (Model 2) 


No. of 
station 


Station completion time 
[ min ] 


TRIA(111.5,144.9,200.8) 


Station name 


Edge cutting station 


Edge cleaning and 


: TRIA(119.2,154.9,214.6) 
sequencing 


Panel production TRIA(368.2,478.6,49 1.5) 


Panel cutting TRIA(226.1,295.5,409.1) 


Profile spot welding | TRIA(174.05,227.5,233.3) 


Profile tig welding | TRIA(212.7,279.3,386.7) 


Section spot welding 


TRIA(284,367.4,516.9) 


Section tig welding TRIA(278,361.4,506) 


Grinding TRIA(85,111.3,154.7) 

The transportation and failure times given in Tab. 6 and 7 
are also valid for Model 2. When the simulation model is run 
for 10- day period under the assumption of two shifts, the panel 
line has produced 16 major sub-assemblies. 


Comparison of Model 1 with Model 2 (Step 9) 


Tab. 11 shows the comparison of Model 1 with Model 
2. In both the models numbers of replication, replication 
lengths and working hours per day are the same, whereas the 
numbers of major sub-assemblies are different. This difference 
demonstrates the effect of the module mounting on the panel 
line throughput. 


Tab. 11. Comparison of the two applied models 


Model 1 | Model 2 
Number of replication 5 5 
Replication length [days ] 10 10 
Hour per day [ hours | 16 16 


Number of major sub- assemblies 11 16 


CONCLUSION 


- In this study the simulation model has been created by 
determining the processes performed on the panel line. The 
required data achieved from the process analysis have been 
put in the model. The simulation model has been run for 
10-day period and as a result the panel line has produced 
11 major sub-assemblies. In the next step, it was assumed 
that the matrix module structure is assembled on flat panel 
assembly. In this case, the completion time of section spot 
welding has changed from 501 min to 281 min. When the 
simulation model has been run in this case, the panel line 
has produced 16 major sub assemblies. This way the panel 
line produces 5 blocks more during 10 days, that means 
that its productivity increased by 45%. 

- Therefore the application of assembling matrix module 
structure increases the panel line productivity by 45%. This 
is a good result in terms of shortening the ship production 
cycle. 
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The Ship Handling Research and Training Centre at Ilawa is owned by the Foundation for Safety of Navigation 
and Environment Protection, which is a joint venture between the Gdynia Maritime University, the Gdansk University 
of Technology and the City of Ilawa. 


Two main fields of activity of the Foundation are: 


> Training on ship handling. Since 1980 more than 2500 ship masters and pilots from 35 countries were trained at 
Itawa Centre. The Foundation for Safety of Navigation and Environment Protection, being non-profit organisation 
is reinvesting all spare funds in new facilities and each year to the existing facilities new models and new training 
areas were added. Existing training models each year are also modernised, that's why at present the Centre represents 
a modern facility perfectly capable to perform training on ship handling of shipmasters, pilots and tug masters. 


Research on ship's manoeuvrability. Many experimental and theoretical research programmes covering different 
problems of manoeuvrability (including human effect, harbour and waterway design) are successfully realised at 
the Centre. 


The Foundation possesses ISO 9001 quality certificate. 
Why training on ship handling? 


The safe handling of ships depends on many factors - on ship's manoeuvring characteristics, human factor (operator 
experience and skill, his behaviour in stressed situation, etc.), actual environmental conditions, and degree of water 
area restriction. 


Results of analysis of CRG (collisions, rammings and groundings) casualties show that in one third of all the 
human error is involved, and the same amount of CRG casualties is attributed to the poor controllability of ships. 
Training on ship handling is largely recommended by IMO as one of the most effective method for improving the 
safety at sea. The goal of the above training is to gain theoretical and practical knowledge on ship handling in a wide 
number of different situations met in practice at sea. 


For further information please contact: 
The Foundation for Safety of Navigation and Environment Protection 


Head office: Ship Handling Centre: 
36, Chrzanowskiego Street 14-200 ILAWA-KAMIONKA, POLAND 
80-278 GDANSK, POLAND tel./fax: +48 (0) 89 648 74 90 
tel./fax: +48 (0) 58 341 59 19 e-mail: office@ilawashiphandling.com.pl 

e-mail: office@portilawa.com 
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